{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simplicial Complex in Three Dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We dsecribe the data structure of the simplicial complex associated to a three dimensional trianglulation give by `node,elem` . The `node` records\n",
    "the coordinates of vertices and `elem` is the pointer from local to\n",
    "global incices of vertices. See [Basic mesh data structure](meshbasicdoc.html).\n",
    "\n",
    "A brief summary.\n",
    "- `edge`: asecond ordering, i.e. `edge(:,1)<edge(:,2)`\n",
    "\n",
    "- `face`: asecond ordering, i.e. `face(:,1)<face(:,2)<face(:,3)`\n",
    "\n",
    "- `elem`: either the positive ordering or the ascend ordering. The default one is the positive ordering and the asecond ordering is mainly used for edge and face elements. \n",
    "\n",
    "- Use `[elem,bdFlag] = sortelem3(elem,bdFlag)` to change the ordering to the ascend ordering. Note that `bdFlag` should be switched together. \n",
    "\n",
    "> The multigrid solvers use the original ordering of `elem` obtained from either uniform refinement or bisection methods. So let `elemold=elem` before sort.\n",
    "\n",
    "- Examples on the usage: `Poisson3RT0; Maxwell; Maxwell2;`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic data structure of a mesh consists of `node` and `elem`. The corresponding simplicial complex consists of vertices, edges, faces, and tetrahedron. We shall discuss three issues \n",
    "\n",
    "- *Indexing* of simplexes\n",
    "- *Ordering* of vertices\n",
    "- *Orientation* of simplexes\n",
    "\n",
    "The indexing and ordering are related and the ordering and orientation\n",
    "are mixed together. However the indexing has nothing to do with the\n",
    "orientation. The indexing and ordering are the combinarotry structure,\n",
    "i.e. only `elem` is needed, while the orientation also depends on `node`,\n",
    "the geometry emembdding of vertices.\n",
    "\n",
    "For indexing, ordering and orientation, there are always local and global versions. The relation between the local and global version is the most complicated issue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Indexing of Simplexes\n",
    "\n",
    "The indexing refers to the numbering of simplexes, e.g., which face is\n",
    "numbered as the first one. There are two types of the indexing: local and\n",
    "global. Each simplex in the simplicial complex has a unique index which\n",
    "is called the global index. In one tetrahedra, the four vertices and four\n",
    "faces have their local index from 1:4. \n",
    "\n",
    "In the assembling procedure of finite element methods, an element-wise\n",
    "matrix using the local indexing is first computed and then assembled to get a\n",
    "big matrix using the global indexing. Thus the pointer from the local\n",
    "indexing to the global indexing is indispensible. For bases independent of\n",
    "the ordering and orientation, e.g., `P1` and `P2` elements, this pointer\n",
    "is sufficient, otherwise, the inconsistency of the local ordering/orientation\n",
    "and the global ordering/orientation should be taken into account.\n",
    "\n",
    "### Local indexing\n",
    "\n",
    "The tetrahedron consists of four vertices indexed as [1 2 3 4]. Each\n",
    "tetrahedron contains four faces and six edges. They can be indexed as\n",
    "\n",
    "    locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];\n",
    "    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];\n",
    "\n",
    "In `locFace`, the i-th face is opposite to the i-th vertices and thus\n",
    "this is called _opposite indexing_. In `locEdge`, it is the\n",
    "_lexicographic indexing_ which is induced from the lexicographic ordering\n",
    "of the six edges. The ordering of vertices of each face or edge will not\n",
    "change the indexing. For example, the following `locFacec` and `locEdged`\n",
    "has the same indexing as `locFace` and `locEdge` but a different ordering\n",
    "of vertices.\n",
    "\n",
    "    locFacec = [2 3 4; 1 4 3; 1 2 4; 1 3 2];\n",
    "    locEdge = [2 1; 3 1; 4 1; 3 2; 4 2; 4 3];\n",
    "\n",
    "Indeed any permuation of each simplex will represent the same simplex and\n",
    "will not change the indexing. The ordering of vertices will affect the\n",
    "orientation and will be discussed later.\n",
    "\n",
    "For a face consists of three vertices [1 2 3], there are two indexing\n",
    "schemes of its three edges.\n",
    "\n",
    "- Oppoiste indexing:        `[2 3; 3 1; 1 2]`\n",
    "- Lexicographic indexing:   `[1 2; 1 3; 2 3]`\n",
    "\n",
    "Each indexing scheme has its advantange and disadavantange and which one\n",
    "to chose depends on the consideration of ordering and orientation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global indexing and vertex pointers\n",
    "\n",
    "Each simplex in the simplicial complex has a unqiue index. It is\n",
    "represented by vertices pointer from the local index to the globa index\n",
    "of vertices.\n",
    "\n",
    "The matrix `elem` is the pointer from local to global indices of vertices\n",
    "of tetrahedron, e.g. `elem(t,1)=25` means the first vertex of the\n",
    "tetrahedron t is the 25-th vertex.\n",
    "\n",
    "Similarly the `NE x 2` matrix `edge` records all edges and the `NF x 3` by 3\n",
    "matrix `face` records all faces of the triangulation. These are vertices\n",
    "pointers. We shall discuss the elementwise pointer from the local indices to\n",
    "the global indices for edges and faces.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem =\n",
      "\n",
      "     1     2     3     7\n",
      "     1     4     3     7\n",
      "     1     5     6     7\n",
      "     1     5     8     7\n",
      "     1     2     6     7\n",
      "     1     4     8     7\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4wEFCgMgS27s0wAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNS1KYW4tMjAxOSAwMjowMzozMvhyhy0AACAA\nSURBVHic7d1rdFXVvffxf0IMBrlqDCo3RW5HpfIobby0KFocFssotYoNlpa0HGqPbT0YfMLQ51iw\nDg+xRbyUdhSlqVIfKtojFo3VoHirFgWfKGpNBCwiYIIa5E7I5XmxwmKz987ea629LnPO9f0MXyR7\nb3YmDpJffnPOtWZeR0eHAAAQtfyoBwAAgAiBBABQBIEEAFACgQQAUAKBBABQAoEEAFACgQQAUAKB\nBABQAoEEAFACgQQAUAKBBABQAoEEAFACgQQAUAKBBABQAoEEAFACgQQAUAKBBABQAoEEqK6lpWX7\n9u1RjwIIXEHUAwCQxapVq8rLyy+55JIePXoUFRUlPnXjjTcOHTo0qoEB/iKQANV99tlnTU1Nf/7z\nn5MeHz169N133x3JkIAgEEiA6tasWfOnP/1p4MCB+flH5thnz569dOnSggK+hWEO/jUDqtu/f/+1\n116b+MiSJUvKy8tHjBgR1ZCAIBBIgOpuueWWxE83bNiwatWqZcuWRTUeICB5HR0dUY8BgAuXX375\n0qVLTzzxxKgHAviMbd+AThYvXlxaWkoawUg0JEBReXl59sfW92l7e/uQIUNqa2tHjRoV3biAoNCQ\nAOXk5eUlppH9yLJly3bs2EEawVQEEqCNRx55hMtgYTACCVBLUjdKtHLlymHDhoU5GCBMBBKgk8LC\nwqiHAASFQAJ0cswxx0Q9BCAoBBKgqLE3jBWRkjElI68aaT/IvYJgMAIJUELeYdUN1WNvGDv2hrG7\nt+4WkcEXDZbD4SQi8+fPj3KUQJD4bQuImLWLoWZLTeP+RvvB3Vt39xrQK/FlY28Ye9vs2wYOHMi1\ngzAVDQmIRmIlqm6oTkyj66+/PvGVvQb02r119/XXX9+4v7G6oTr1KiXADDQkIGxpK1Gi+TfPn/mL\nmRWjK0rGlIhI2YyyxfMW289WN1Rbb0JVgmFoSEBIMlSiRFYaJT048xcz59981OoRVQnmoSEBwbIz\nw2o2PqIqwTA0JCAoVoOp2VJjVSInf8SuRxWjK+a+OTfxqdSSZKEqwRg0JMBnSZWoq6m5VGkn65yw\nq5Icvi84oCMaEuCbxCWi3Cfolt2x7D/+8B+Jj3RVkizVDdU1W2qoStAXDQnIlS+rRKn1qHlDs9s3\nsfeFC1UJGqIhAR4l7ZrzMY1SF5BsmUuShaoETdGQANeyXkgUOaoSdERDApxyeCGRW273MjgpSRaq\nEvRCQwKyCO5CIslhZ51DVCVohEACuhTJ1FyGBSSbVZLm3DHH4XtWN1T3L+rPJbRQHIEEJPN8IZFb\nQdejRFQlqI9AAo4IsxLlmEZuS5KFuw1BZQQSEF4lsmVNo9/+8LdJV8X6iKoENRFIiK9Adyt442QB\nyeatJFmoSlAQgYQ4ivZCojCXjjKjKkEpBBJiJPypuVT+plEuJclCVYI6CCTEgvr3VogWVQkqIJBg\nMhUqUaLM9cjVApIt95Jk4QwLRI5AgoEU3K0gKi0dZcAltIgQgQSjMDWXOy6hRVS4uSpMENBtT32U\ntR4lztelHs2XmfPbrTrHjVkRPhoS9KZFJXI7WefhaL4gUJUQMhoStKR+JQpZECXJQlVCaGhI0Ima\nuxUy02IvQ2ZUJYSDQIIetJiaS+Uwjbxt+E7i1/7vrrABD0EjkKA01S4kijmqEgJFIEFRmlaiROFP\n1gVdkizcbQgBIZCgFmMqkQFLR5lRleA7AglK0HG3gl98WUCyhVOSLFQl+ItAQsQMmJpLlUs9CvRo\nviBQleAXAgnRMGZqLlXkk3VhliQLN2aFLwgkhM36sRXDqTnjsS8cOeJODQhJ0r0Voh5OUNzWI38X\nkGzB3bghM3tfOHd2gAc0JAQrVrsVIp+sUwRVCd4QSAiKkbsVNBL+SlIiLqGFBwQSfGbwboXMqEep\nqEpwhUCCb+K8W8FbGgW0gGSLtiRZqEpwjkBCrmK1ShQ0t0fz6YKqBCcIJHhEDtl8nKzz/Wg+FUqS\nhaqErAgkuMZuhUSe0yjo+To1UZWQAdchwSkOadVUVNckdYVrldAVGhKyYGouA3bWecaNWZGKQEKX\nmJrLTKM0UmclKQmrSkhEICFZbC8kciXHNIrnAlJaVCXYCCQcQSUymLIlyUJVghBIECqRexpN1mmE\nMyxAIMUXuxW8CS6NtDuaLwjsC48zAklpe/bsuf3225ubm1taWvbu3fud73znmmuuyf1tmZqLVlQL\nSIrP2tm4hDa2CCR1tbe3T58+/Z577hkwYICItLa2Tp06dcuWLbNnz/b2hkzN5Y7JutBQlWKIC2PV\n9Ze//OXcc8+10khECgoKFi1adNttt3l4K+sixJotNWYfjhc0A9JItYtkM+MS2rihIanrxRdfbGpq\nSnzkhBNO2Ldvn/N3oBLBAFSl+KAhqev8889/9NFHf/7znx84cMB6ZNmyZU7WkJLu8UMl8oVf9Sjy\nK5D0KkkWqlJM0JDUdfXVV8+bN+++++5buXLl/fffLyKPPvroQw89lOGPsFshIAZM1hmAqmQ8Akld\nhYWFL7/88ve+971Vq1ZNmDBh6NChr7/+eu/evVNfydQcXNFlu10qNuCZjSk75eQlaGxs7NWr169+\n9at+/fpt2rTpjDPO+Pvf/576YnYrBCqcemTq0XxBqG6ortlSw/SdeQgkhaROkZ999tmPP/747Nmz\n6+vrr7766qampkmTJn322WecBBEaf9MowwKS70fzZabjSlIiVpWMxJSdHk488cTly5fn5eU1NzcX\nFxdThgDhxqzGoSGpIsMvenYfqtlSc+6Ec8ddPS7MgcWZ2XsZdC9JNqqSMWhIerB3K/To0+PUs06N\nejix4HsaRb7h22BUJTPQkHSyZ+eebRu2XfCtC6IeCAxhTEmyUJV0R0NSXWHvwqLiondeeWfHlh1v\nrnpzxp0zuvfoHvWgzGf2ZJ3BOMNCazQkVaT95ikZU9Kyq6X81vKDrQcHnjmwYknFSaeeFP7Y4iZW\naWRYSbKwL1xTNCTVlYwpuXv63cKpRZpjASlkXEKrIwJJIda3TeKvdU11TXYOlY8oJ5NCEEk9ivZo\nPn1v3JAVdxvSC1N2yuno6LC+eWq21CQ+Xt1QXT6ivHxEeUTjioVYTdbFBJfQaoRAUlrNlprEBLLu\ny0AmwXdGriQlYlVJCwSSfqhKAQmuHrGApAKqkvoIJNUllSQLVcl3TNYZX5IsPlalXbt2zZ49e/Lk\nyeXl5eXl5evWrcv9PWOOQNIYVQnwwJeqtGPHjvPOO2/YsGErVqyorq4eN27clVde6eMg44lA0kDa\nkmShKvmCemSJSUmyWN87njNp6tSppaWl1113nfXprbfeunXr1tbWVv8GGEds+zaBnUnsC/cg6DRi\nAUll3q5VevbZZ1etWvXqq6/ajzzzzDN79+4tKOAnak7436eHmi01EwdNzJA31lPEklsqdCOljuYz\n+Jqkrni4MWt1dXW3bt1KS0vtR84444xABhczBJJR7FgikzQS8tF8SMtVVVq5cmXv3r03bdq0ePHi\nL774Yvv27d/61rd+9KMfBT9MwxFI2shakmzM4DmkQj1SUAxLksXhjVkPHDiwd+/ewsLCpUuX3nnn\nndYjZ5555qefflpZWRnaaI3EpgYzsdkhq3DSiAUk7XS1L9w+J7OoqEhEmpuby8rKrKeOPfbYqVOn\n3nLLLZs3bw57uGYhkHSSYbtdWuwLBzxI2hfe1QbxUaNG2R+PGDGira1t5cqV4Y3SRASS4ahKaTFZ\nl1ms9n93xdUltIWFhSLy1ltvBTwowxFImnFbkixUpUShpRHzdbpr3N+Y4dnErGpvbxeRkpKSwMdk\nNAIpLqhKcIWS5EpjY6OIXH755VEPRG8Ekn68lSQLVYnJOgThiSeeuPjii7/2ta9FPRC9EUixE+eq\npGAaRXs0X2aUJBEpGVNSMqbLibiHH3742WefvfHGGw8cOLB8+fIwB2YkAklLuZQkC1UpaCwgmaGp\nrqmpriltLH366acFBQVbt26dNGnSa6+9duKJJ0YyQpNwYWx8xe22DgrWI/XF9iJZ+3e16obqyimV\nItJU15T0mhNOOOGaa64Je2RGI5B05fzGDZnF5LYOpBEcSv12qFpelTSX4OpOrHCOQAI3ZkUmMSlJ\niZUo7QtqttSIyMRBE0mj4BBIGqvZUjN90vSq5VW+vJvBM3jh1yMWkHSRNYes11hphKARSHpLndfO\nkXkzeEzW5c7IkuTw3zlpFCYCCckMrkqAk0qEqBBIevNra0MqM6qS4vVIqaP5MjOgJHn490w9ChmB\nhC7pXpWiSiPnC0gczRcCz5WINAofgaS94EqSxYyqhNzpVZKYmtMRgYTsdKxKik/WITi+/P5EPYoE\ngWSCoEuSRaOqRBoFROWS5GMlIo2iQiDBBR2rUsi4AilkTM2ZhEAyRDglyaJ4VaIeBUqdkhTQP0Lq\nUYQIJHihbFUijYwXaCUijaJFIJkjzJJkUbwqhS8m83VRlST+sRmPQEJOlKpKetUjlY/mU0poq0TU\no8gRSEYJvyRZVKhKeqWR7kIoSSHvViCNVEAgwR+cYQG/8K8otggk00RVkixRzeCpUI9isoBk870k\nRbiBm3qkCAIJ/gt5Bk+FNEIuoq1EpJE6CCQDRVuSLKFVJdJIX1zTiiQEEgKkwmYHBMfbrJ1SOUQ9\nUgqBZCYVSpIl0KqkTj2K2wKSN6r9dkIaqYZAQhiCqErqpJEHGh3Nl5mTkqRUJYLKCCRjqVOSLEpd\nQhu5mBzNp1olSkQ9UhCBhFD5VZW0rkeGSS1J6lci0khNBJLJVCtJltyrkmppxAKSRf0cguIIJESD\nDXgmqX+svvwxbSZjqUfKIpAMp2ZJsnirSqrVozirGF1hfbBg/YLF8xZHOxiHSCOVEUiImKuqpGAa\nxXO+zoqiBesX2I+oc3BfBpVTKkkjlRFI5lO5JFnYgKeLxEoU7Ui8aaprinoIyIRAgiqyViUF61FM\nOMwhxUsSk3XqI5BioWZLzfRJ06uWV0U9kCwynGFhUhppdDRf6tScpkgjLRBIcaHRZIVGM3imLiB5\nnppTvCRBcQQSFJU4g2dSPVKcMZUoEfVIFwRSXKi/tSGVXZUM+/moIB93K6hWkkgjjRBIUN3Iq0Ya\n+Wu7CnTfNQfDEEgxomNJ6pys+4WISMXoCqV+bmq9gBRoxqtTkqhHeiGQoK6kpaMF6xc4/zFaNS3N\nlsLKpZU+Dk9HsapEpJF2CKR40bEkJbJ+jGaoSok5tHr1UU+NH3/k2RgmU8jTnuqUJGiEQIKiMuys\n66oqVU2rSgqhRIlPjR9fFVUmhXw0X6wqUSLqkY4IpNjRoiRl3eedWpUyp1GS1atzzSTPC0jhHM2n\nQg5FWJJII00RSNCYXZVKxpQ4TyNL7pmkJnYkQl8EUhwpXpJcXQa7YP2CqmlVTXVNIiVuv5ABmWSv\nitl34igZU6LO3yiSkkQ90heBBLW4vSnD4Zm6kjP7NYnIu83uYknTTErNIfsvHvO9G6SR1vKjHgCi\nUbOlxj5wWl+J60bvNpe829wZS66sXp1+j3gG0V6BZP2tm+qamuqarL91YgyvXn3kP7d/ryBYJSnq\nUUAPNCQoJPd71tmZ1FVVamqUlStFREaOkq9+NZcvFQ1riejMfo66oKb9zzPqke5oSPGlWkny6w6q\ndlVKbUttbfLLX8qfHpY/PSxr1uT+pUJVMbqiYnRFaiXKTIWeFE5JIo0MQCBBCR7SKPPP2bQzeEse\nkHfe9fJuiUKer1t7z1orikrGuMihRCpkEuAEgRRrqpUkt7Ju9U6sSq+vkT8/4vF9/OL8aL6558y1\n/ht7w9gF6xd42NeeyPhMoh6ZgTUkRC/Q446sVnFmv6bhF5V0dMipQ+RfmwP6Uv6Ye85cEbFL2JMP\nPunqmt+uRLueFOj+b9LIGARS3EV+TVI4h+9N+H5J7UNNInLvPSU/vyHor+aFlUOSEEUe7N8nj6+Q\nhgbZ2SynDJCLL5KvlPoyOiAMBBLMt2SJvPW2DBhbct9v5MrzmkrGeFmJsVSMrhj3/XEbn9tYUFhw\ncO/B+lfrB501aMyUMbkML6kSebb5X1JZKY2HV83eelueflrGfU3m3ZbjG/sjoJJEPTIJgYQoS1II\n9WjdWvm/D4uI/Pg6GTlSVq4tmTS2SUTkmx5j6aWHXrI+KCgsuOwnl3lOI4eV6KMXP3L4hnfe2ZlG\nhYUyeLBs3CgdHfLSy7J8uUyZ0vmaqmlG7QInjQxDICEyuafR+PHZ9yMsWSLtHdItX5qbZckD0rxT\nrIb0h9ub+vYpuegiGTZcxo938UWvnX/tgb0Hjj3u2BEXjOjo2eF2zB6m5pwsIG3cKO/9U0TktNNk\n4ULp00dqn5U7/ltE5PU1nYG0erW7v6nvOJMCmRFIEFFgJcmDyqWVTnaOtbaJiLS1y8MPH/V4yZiS\nu25quuvwrgfnvWH4ZcOtDzrEXRr5NTWX1gcNnR9cdJH06SMiCatHeUF8wehRj8xDICEa4exlEJGe\nPaVP7yOftrXLnj0iInl5MnxcyezZcmY/p6tKFaMrRl862u0AEivRb3/4W7d/3KHLvyGXXCoiUlgo\nX+yUjz8+ssf9y18O6Gt64VdJIo2MRCChU5glKbQ0EpG77jrq000b5UczRESuuEIqKkQOX6uU4RTa\nJDvW7/jgtQ+OO/647Q3bx31vXI/BPbp6ZVIlCvpovsLCzg9u+T/y7uHrf8ePlyuvDO5rRqNySiVp\nZCQCCRqrXFo5fnyu1+iMH5/9ZHTbR+98tHP7zgtmXiAiQ+qH3F1298zfzyw+qzjxNV2tEoVzNJ+I\n5CfM0b3xhrywWiZcFs5XdiT3kmTf4xyGIZBwRM2WmumTplctD/aSfn/rkedMys8XERk//sjqUdqT\n0e2jVy3ld5f3HdnX+rjvyL7DS4cvu3nZz/76M/HpQqK0nOzdSHTvfbJnj9TWyn33yZ49cu99cskl\n0q0g4h0NfmGyzmAEEo4S9O+eQUzWucqkoacf+eGemEYWuyp19cfvLrs7MW96Ffd694V3A92tMPii\nwU52fv/mN/LOehGRX/9aevaSnj3l29+WlX+VD/8le/bI22/L/zpHRJkTkjyXJNLIbAQSTOChJ6Wm\nkS1rLFnsSnTt/GvtrXdRaWmR+gYRkdffkEsuERHZu1c++aTz2eN6RjYwwDluroqjBHq71UD3MlQu\nrXQ+JZUhjSwZ0si+7encN+d+5dtfEZH8btF/H40+vPvv3nvk6afl+edk7lzZf0BE5LgecvrpEQ4t\nPQ9nUlCPjEdDQkhC2Fln9ST706TClBhXOc5c2bNze5v35ufnn/7lYH/eO+l/EyZIzVNS95Z8sUvu\nvPOop266Sbp1y57BiiON4oBAQrIQ9n8vnrc49UFf4sr+mVs1rSq1MPn7EzlvT977r7z/9Zlfl97Z\nX5wjJ5l0xx3y4EPy+OPS0tL5yOlDZca/y3nnKZpG3LgBSQgkhMGqR4k5tHDhUS+YNetISvmbTEHY\n9sa2tkNtz//h+Yu+f5G1BTwEWTOpqIdcd51cd51s3Sq7dsmAAdK7t4iD+Un1UY9iIq+jw/XNuBCC\nvLy8aL8DfSxJdholhVBXZs3yJ5Ny1NUy0uy/zN5Wvy2/W/5p55zW7YRuDt/N+dF8SZ588MmyGWX2\np27PRlI/jRbPW5y5JCmVRhMHTeRnZnCiX4xFHDhPIxFZuDD9nJ4iep7Wc8TlI4ZNGOY8jXzk794N\nQCkEEtLza7udtZPKeRpZVMikBesXJN24Ye6bcwO62MgVK5Ps/5IkPqVFGmXebqdUPULQWENCgLyl\nkWXhQpk1a3Hkc3f2NUkqRJEtzL0bESKN4oZAQpd82W7nLY3sP6tCJqnMjOxhux0sTNkhKPNvnp9L\nGllUmLtDJKhHMUQgIRPPK0luL8L/oEFqnpJFi+SPf5SXX5bWVg9fExpLWkkijeKJKTtE74kn5IUX\njnz61lvy2qvy8xvk2GMjGxKA8BFIyMLDSpKrevTcc51p1Ke3jPo3qa+XnTtl+yfy17/KlCmdr1k8\nT++VpKCP5ouQjzfdsFeSqEexRSDBZ64ug21rlRdfEBHp0UP+8z+lbz/ZtUtuv10OHZL16zsDaeFC\nmTUr0CEHLrSj+UIT3E03SKM4I5CQXXB3t3vnHdm9R0TkvPOkbz85cEC6dZP/+i9pbe08QA8Kyvzb\nRuJTbJKEKwQS/OT2lt5bt3V+UNRDquZLY5N0dEhxsUyaJF/6UiAjRI7c3nTDVSbVP1ZfMqbE48ig\nP34LhSNOttt5OGCi+fBU1lNPyieNYt0k7NNPpbpa/rHGyzgRKFdpZHG+cb9idEWg95iH+ggk+MPb\ncUcH9nd+kJ8v375SZs+WCy/sfKTmSWlj87dKPKSRxdXFZFXLq6ZPmu7ly0B/BBKccntNkpOdCPbR\n2mPHyrivyYABctVVUnyCiMjuPbLjU6fvg6B5TiNL1kyiHkEIJPgitR45bEu9Dx9td9JJRx7sf/jj\nffvcvRsCkmMaWTJkUlIaUZJii00NcCHtdrtcziYfMVxqa0VE/vWvzkfa22Tb1s6P+xuxvG3wRUhJ\nWlpkVW2axwefKmedGfpooCECCVEaNlyO7yefN8vbb8vKlTJypLz2mjTvFBEZOPDIhJ7WzLsIqSvb\ntkntqjSPl5ZmCqS0k3VWSfrjyj/6OT4ojyk7uJO0kpShHs38xUwnyz9TrpGCAhGR55+X3/1O6upE\nRPLz5aqrRJQ5PTbOnO9H2L6984MePeS4hP+OLery3Vg6QiIaErzLOlk38xczZ83KsvwwcqT8+wx5\nqkY+/lja20VETj5ZvnuNDB5CGqnC4QLStu0iIvl5MneuHHNM+vdxvkWFkhRDBBJcc3XjBieZNGKk\njBgpBw9K4ydyQrEcd5wI3UhD27eKiBSfKPl5sm2btLXJKSdLt65/xlCPkIQpO3jkfC+Dw7m77t1l\n8BDSSGPbPxERaWmR234pv/qV3HWXzJkjTz+d/sVO0ojtdnFDIMGLkjElrtbqHWaShTTS0c6dndv0\nd+6UXbs6H2xtk2eflWeeSX5x1bQquhFSMWUHj5rqmly93pq7sz9NvUV04itzGhmiYO9oGDhQysqk\nTx95+SV55lkRkdXPy8UXSfeE062c/+NhJSlWCCS4Nn3S9KrlVR5OkrWTZvG8xamFiRxS1qxZ2fc1\njBwp8/9bROSYws47tV/+Dfnn+/LRR3KwRRobO3epCEtH6BqBBHesNBKR6obq8hHlC9Yv8PAm8cke\nA66KtU63yvqyXbtk/34RkeLiI0eHFBfLRx+JSOchIyJS/1i92zSiJMUHgQQEKD5XxT7xROc1ZN/8\nplx6qYhIa6ts2tT57CmnRDYwaIRNDXDBrkeW6obqitEVEY4H6hg6tPOD55+XVavkn/+UJQ/Izp0i\nIqecLP36iXiqRxa228UEgQSnktII8eFkk+R558mgQSIi+/bJU0/J4sXyfr2ISEGBXPs9EZGJgzym\nEeKDKTvkJJeVpAgte2BZ6oNlM8rCH4kTK1asmDBhwqZNmw4dOtTR0dHS0jJo8KCBAwaGOYasFzgf\nc4z8eKb87Rn5xz+k9fBBVkOHytVXy0kn+XCGCCtJcUAgwREz6lFiDl1//VFPLVp05FkFk6m2tlZE\nJk+eLCLbtm97fc3re3btGfVvo8IcQ9ZMOq6nfOc7Mvlb8tnnsm+f9O8vRUUiIrNmeZ+sQ6wQSMgu\ncxrpUpKWPbAsKYQSJT61aNEyBTPJdsrJp4hIwwcNIQeSOLsRVLcCKUk4N8THNKIkGY81JMRC5jRK\ncv316ef0IrdixQr74zzJi2QM3HQDwaEhIQsnk3WKlyRXaWS5/nofelJAFyF9/vnnIjJg4ADf39kh\nVzfd8PcyWEqS2QgkZGLA0pGHNLLknkm5X4T03FPPpT5YX1/ft2/fs846K+sfD27vhsObbnBTBrhC\nIMEfapYkz2lk8aUneZCYQ0XFRfs/3V9U3HnI3QsvvNCjR49x48bl53c53x7m3o3wZ+QoSQYjkNAl\n3etRjmlksTLJj+E4YkXRpVdcan365INPFhUXWZkkhyPkrbffeumll8aeO7ZnrzRnvKuzd4N6BLfY\n1ID0PKSRvjduaG+X996T55+X//kfqX1W3nm38+za0Dz31HPWf5decamdRomsWLLKzdlfOnvf3n0v\nv/LyoZZDSS9TZ+9GoGnEjRtMRUNC3O3fJ0//7cjpCdtFGj6Qt9+Sb18pRcdm/JN+SKpEXbGuQJLD\n03Gnf+n0bdu2bd6yedjpw+zXRLV3A/ALDQlpeJ6sU6ckOf/d/9XXOtMoP196Hp4Da26W1Wn2E/gm\nayUSkScffLJsRtnkyZN79+799N+e3rt3r4iUzSgrm1G28e2N+z/dv3v3bvvFuezd8L0nMVkHbwgk\nJNN96cjm5Af0nj2dt1wrLJRrvis/+IFMvKLzqc0fdd4CJ/eFKFtSDmUtRiLS3t6+a9eugwcOHjx4\n0H7w1DNPLSouen/t+9anue/d8DGTwkkjZu2MxJQdfKbmdruuNDeLdIiIDBsux/cTETntVOnTW77Y\nJe3tsmuXHH+8P1/I4dRcqvz8/OOPP37w4MHHHx7KoZZDO3bs6Nmz56RJk6wg8WvvBnN3iBaBhKMY\nU48camuX4mIRkQGHLzM9dEj2HxARkTzp3cv7O1tXxdobuD1Eka20tLSurq6lpaVv374d0rFp46be\nvXuXfqW0W7duZTPKPJSbdW/KoRYZMKDz/tw+CnOyjv3f5iGQcIRfaaRRSTp1iJw65Minra1Su0pa\nWkREBg6QgmO8v3PzhmZras7bH7cWkKyPu3fvXlpa+sUXX+zavUs6ZOSokcf3817cGhvlH6+JiLS2\n+RxILB0hRwQSOsWtG6XasUNWPSeffyYiUlQkF1/s5U0Sr2nNpRWl6tOnT58+fZIedFWP2ttl80ey\nenWXL1j2gGazdpQkwxBICIQKJWnRIheLK2/+P1nzj87Lj4qL5RvfkN69j7xP1D0BwQAAD5RJREFU\nVqlTc/V/qXc1Ws8c/h1XrZING6StLdP7OPmbdoV6hNwRSBAxsR45X1xpb5dna2XjBhGRvDwZO1bG\njpWu78uTzPNuhZAdPJApjXIUYRpRkkxCICGoNFKhJDmx+oXONOreXa64Qk4+2dGf8mW3QlcSF5D8\n8qWzxbqIdvsn8t67/r434A8CCbH2SaO833kxjxSfIB9+KB9+eOTZ0aOlV8pGO10qURJ7/0Jens+B\nFPlkHSXJGARS3AU6WVfdUF05pbJyaWVA759Z2YyyRYuyXDH60ebO65BEZOs22brtqGcHDpKHHpKy\nGWVr71kbaCXSV+RpBJNwp4ZYC2HpqKmuKdD3z6xsRlnmhfqdOzM9u/KvIod3sjm/t4IEdjRfWrns\nRAjifSLBjRvMQEOC4TL3pMsuk8suS/+U9QPaWstZe89aV180l6P5XC0gebswNsO7uXo99Qj+oiHF\nVzg761S43WrWnpRq0aLO25gGMyITqJZGlCQDEEgxZd4+78xcZZKVRkEOR3uqpRHMwJQdAqfI/m9r\n7s7+NPVs78RXhjWoZM899Zzbr+5k70ZWZmQw2+10RyDFUdzqkc3+mbvsgWWphUmFn8jWUeVuucqk\nvJRpEbdpRD1CQAik2IkkjRQpSRZrF4AK8eMj55k0YriMGH7kUyuYnd/FTvE0oiRpjUBCXNi70QyL\nIpuHubvEbqTdnVVhHgIpXiKcrIuqJGmUQ7nfMchVJiXN1Nk7yDOMQfF6ZKEk6YtAipG4LR1FNTUX\n5lWxqXLZu2E90lVV0iKNoDUCCeEJpyRFXolyuSrWFznu3XBSlRRHSdIUgRQXcahHuv8Y9Z3n/xWp\nVYl6hBAQSLGgThoFUZIir0S+COLIiRzZVWntPWu1SyNKko4IJOjKjBxSnLdb+QHeEEjmU6ceWXIv\nSUzNhcmarJt/83wRmXPHnKiH4wIlSTsEkuFUS6NcUInCZy8dWVE0/+b5emUS9EIgIQJuS5LxlUjB\nBaS05twxR8eqBF0QSCbTvR7pWImivQjJX2l31tlVSXSIJWbt9EIgGUvxNMpQknTMIVvkFyH5JfM+\nb2bwEAQCCQoxfmrOMFrM4FGSNEIgmUnxemSxS5LWlSh3Ci4gOb8MlqoEHxFIBtIijWzcZFo1Hm7K\noHhVoiTpgkAyjRZpZP3wksMliUAyAFUJuSOQEB47h/iZZfNwZnmgcrxnnbJViZKkBQLJKMrWoww/\npJQ6TDZ83s4sD4gvd1ClKsEzAskcCqYRlSi2FKxKlCT1EUgIhKsfRtUN1ZVTKiuXVgY8qMBpfVWs\n7wdMUJXgFoFkCEXqkedK1FTXFMBwwub2qlh1NnwHd9yRUlWJkqQ4AskEkadR7lNzxpQkpKIqwSEC\nCTnx8ZdfM0qSjsI5DVaRqkRJUhmBpL1I6hG7FYwR5tnkVCVkRiDpLfw0Cu6X3Ljt/1ZhAalqWlX4\nZ5NHXpUoScoikOAIlchIUU2T6nWGBUJDIGkshHoUcg7FrSRFK8zJurQinMGjJKmJQNJV0GnEb69u\n6XURUuRpZIt8Bg/qIJBwlMin5vQtSa4uQlJhAUkdkVQlSpKCCCQtBVGP+C01PtSpR4moSiCQ9ONv\nGkVeiVLpW5K0oGYaWUKuSpQk1RBIMaVgDsWHakdOqIaqFFsEkmZyr0dafKubXZIiPHJC5XqUKLSq\nRElSCoGkk1zSiEoEXdLIRlWKGwLJcPrmkNklCQ6FUJUoSeogkLThth7xq6WaolpA0q4eJaIqxQSB\npAfnaaRvJUqlUUlyflVsJAtIWqeRJdCqRElSBIFkDn6FjJDbo/ngDVXJbASSBjLXI5MqUSqNSpKy\nDKhHiTjDwmAEknLy8vKsD2bNmrVw4UI7jXY07njxby9e9YOrrGfNziFThX/HIMPSyOZ7VWLWTgUE\nkkLsKLJNnzTd/njln1f26NVD4jc1R0lCWpxhYR4CSXX9hvUTkb8///emxiZpZKYCLphajxL5OINH\nSYocgaSK1HrUvKHZSiPrd8BhI4e1trV+94ffjWBwUaMkeRCHNLKx2cEM+VEPAI7MuWNO927dC/ML\nox4IvOPIiUDNuWOOHUueWSXJpxHBNQJJXUnHS+f3ze/Wr1tUg4lcdUN1xeiKqEeRnoJH88WqHiWy\nMinHWEJUCCSl1T9WLyK9BvSKeiDIRLWLkGKbRpYcqxIlKUIEkio6OjrSPl7/WP3ae9aGPBg1qVyS\noBqqko4IJA10dHSUjyhv39ke9UDgXWgLSDGvR4k8VyVKUlTYZacQqyQlbrerra0VkVWrVtXW1k6Y\nMEFEzjnnnG4F3U45+ZSoBhmt6obqyimVlUsrox6IokijVGzA0wiBpBwrltasWbN06dIVK1bYj3fv\n3r1///63XnrrT373k9gGkqTs9QCy8nCtEtckRYJAUlRpaWlpaWniI4888siIESM2b96cl5f3O/kd\nvwhrJJwjJ6hHmVGV1McakjYOHTrU3NwsIh0dHdaqUvmI8qgHFQEdtzaEcOQEaeSEq1UlVpLCRyBp\noKysbNSoUV988cW6devOOOOMadOmSUIsRT06QDNswFMWU3YaWLZsWVdPdXR0WJsgYvXbsVJ3ElLh\nqljqkVsOb8zKSlLIaEjaoypFK+tVsUEvIJFGnvlytyH4iEAyRNxWlTRaSQp0AalqWhVplKPMM3is\nJIWJQDIHVSmG2ATvC6qSIggk08SnKmlUkgLCZJ2/uqpKlKTQsKnBQPYdH2q21DTub4x6OLEW3B2D\nSKMg+HjcHzwgkIwVhw14Sm23gzFSL6G1SpJ1YGZeXl5HR0dDQ8O8efP279//3nvvDR069Kabbho/\nfnyUgzYCgWQyqpKpqEdBS6xK82+en7SX0vpV7/PPP+/Xr5+IzJkz55JLLrnvvvt++tOfRjJaY7CG\nZL6Ojo6JgyaauqoU7UpSJBchkUahsauS1Y2SYslKIxGZP3/+0KFDb7zxxq1bt4Y/SJMQSLFgb8Dr\nX9Q/6rEYJfNFSJxZbrbEG/N/9atfPXTo0CuvvBLheAxAIMWIqVUpVtvtqEdRsUpSVwoKCkSktbU1\nrOGYiTWkeEk8comfa9ohjaKV4aqvDRs2iMi5554b4nAMRCDFkR1Lxvx0Y7sdAtJrQK+196zN/Jr3\n3nvvpZdemjp16qhRo8IZlakIpPiKw77wCPm+gEQ9Cl/m+W3rFzsR+fGPf3zBBRf8/ve/D2VQJiOQ\nYs2kqmR2SSKNwmTnUOIUd1dmz55dWFj4xBNP9OzZM4zBGY1AAlUJEEnJIVtqLNkvuPfee7dv3/7M\nM89YmxqQI/4nQsSUqqROSfL3yAnqUdCsKErKoSSpzz7yyCPr169/+OGHrU/XrVt34MCBCy+8MKBB\nxgGBhCOoSq5kuCrWxyMnSKPgdFWJnFi9evW6devuv/9++5G6urohQ4b4NrhYIpBwFN3vNhRmScp6\nNB+U5aQSZVBXVzd16tTzzz//yiuvtB5pb29/4403PvzwQ9+GGEsEEtKgKimCeuSvXCpRokmTJn3y\nySePP/544oO9evUqLCzM5W1BICE9fatS5CtJfi0gkUZ+8SuHbFu2bPHlfZCEQEImVCUPAj2zHK7k\nODWHkBFIyELHqlTdUF05pbJyaWXUA/GOepQL3ysRwkEgwRHtqlKG246pjzTyjEqkNQIJTulVlaIq\nSRw5EQkqkRkIJLijUVUKtCQFdzQf9cg5csgwBBJc4wwLCewiJNLIIabmjEQgwSP17zYU+f5v+I5K\nZDYCCTnRaAYvBDkuIFGPutK/qP/EQROFHDIdgYRcqVyVNCpJpFFaVKJYIZDgD6pSLqqmVfH/LRE5\nFE8EEnyjZlXSoiRpfdWUv9itEGcEEnwW26rkeQGJyTqhEkFECCQEQbVLaFUuSTFPI3IIiQgkBMXg\nqhTcVbHxwdQcUhFICJA6VcnfkpR6Vay3IydiWI+oRMiAQELgDK5KNg9HTsQtjahEyIpAQhhUqEoq\nryQZjEoE5wgkhCcOVckh4+sROQQPCCSEKtqqFFBJcruAZHYaMTUHzwgkRMCwqsSZ5UIlgh8IJEQj\nqjMsIl9JMq8eUYngFwIJUVLzbkOZ5XIRkklpRCWC7wgkRC/kGbwcS1LSRUhxO7OcHEJwCCQoQceq\n5Jbu9YipOQSNQIJCQqtK4a8k6ZtGVCKEhkCCWuJQlXRBJULICCSoKISq5EtJcriApFc9ohIhKgQS\nFKXC3YZ8oUsakUOIHIEEpQValSK/JkkRTM1BEQQSVKd1VVK5HlGJoBoCCXoIqCpVN1RXTqmsXFrp\n8PWJV8VmXUBSNo2oRFATgQRtBFSVmuqanL849Wg+jVCJoDgCCZrR5cas6tQjcgi6IJCgH39vzOpt\na0PmIycUSSOm5qAXAgm6ivYSWpWPnKASQVMEEvTmywyev/u/o6pH/Yv6Txw0UcghaItAgvaUuttQ\nJGlEJYIZCCQYIseq5KokKXLkBDkEwxBIMEfQVSnr0Xyh1SN2K8BIBBJM47kqZS1JmS9CCiGNqEQw\nG4EEA4W/qlQ1rSq4r8VuBcQEgQRjeahKTlaS0i4gubrdg3NUIsQKgQSThXNjVt8n68ghxBOBBPO5\nqkpur0nyN43YrYA4I5AQC4qfYUElAoRAQqw4rEoZSlLSAlKO9YgcAhIRSIgXH6tSLmnE1ByQikBC\nHGWtSqklyboq9skHn8zl61KJgAwIJMSU2zMskq6KdVuPqERAVgQSYi3DJbSpJcleQHKeRlQiwDkC\nCfD/FFpyCPCAQAJEuqhKabfbZa5HTM0BnhFIwBFpq1LF6Arrg7IZZVYatTW2/Xr2ryf/bPLIr4y0\nnqISAbkjkICjJG52SFI1rUpE7v3Rvbt279pYt7HtJ21CJQL8k8c3EpBWUiaVjCmx7qBa3VD96opX\n7//f91uP8x0E+IVAAtJI25CS1NbWfv3rXw9hMEBM5Ec9AEA/Dz30UNRDAAxEIAEAlEAgAQCUQCAB\naWRYW2XZFQgIgQQAUALXIQHppV6QRDcCAkUgAZkQQkBomLIDACiBQAJc2717t4i0trZGPRDAKNyp\nAXBh8uTJH3/8cV1dXVtb2zHHHHPhhRcOHTp0yZIlUY8LMAGBBABQAlN2AAAlEEgAACUQSAAAJRBI\nAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAA\nJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQ\nSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJRBIAAAlEEgAACUQSAAAJfx/kKTm5gy8\nbJcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[node,elem] = cubemesh([-1,1,-1,1,-1,1],2);\n",
    "showmesh3(node,elem,[],'FaceAlpha',0.25);\n",
    "findelem3(node,elem);\n",
    "findnode3(node,elem(:));\n",
    "display(elem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate index pointers for edges and faces\n",
    "\n",
    "One can easily collect edges and faces elementwise. The issue is the\n",
    "duplication. For example, each interior face will be counted twice. The\n",
    "`unique` function is applied such that each edge or face has a unique\n",
    "global index.\n",
    "\n",
    "**Edge and Face**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "edge =\n",
      "\n",
      "  19�2 uint32 matrix\n",
      "\n",
      "   1   2\n",
      "   1   3\n",
      "   1   4\n",
      "   1   5\n",
      "   1   6\n",
      "   1   7\n",
      "   1   8\n",
      "   2   3\n",
      "   2   6\n",
      "   2   7\n",
      "   3   4\n",
      "   3   7\n",
      "   4   7\n",
      "   4   8\n",
      "   5   6\n",
      "   5   7\n",
      "   5   8\n",
      "   6   7\n",
      "   7   8\n",
      "\n",
      "\n",
      "face =\n",
      "\n",
      "  18�3 uint32 matrix\n",
      "\n",
      "   1   2   3\n",
      "   1   2   6\n",
      "   1   2   7\n",
      "   1   3   4\n",
      "   1   3   7\n",
      "   1   4   7\n",
      "   1   4   8\n",
      "   1   5   6\n",
      "   1   5   7\n",
      "   1   5   8\n",
      "   1   6   7\n",
      "   1   7   8\n",
      "   2   3   7\n",
      "   2   6   7\n",
      "   3   4   7\n",
      "   4   7   8\n",
      "   5   6   7\n",
      "   5   7   8\n",
      "\n"
     ]
    }
   ],
   "source": [
    "totalEdge = uint32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...\n",
    "                    elem(:,[2 3]); elem(:,[2 4]); elem(:,[3 4])]);\n",
    "sortedTotalEdge = sort(totalEdge,2);\n",
    "[edge, ~, je] = unique(sortedTotalEdge,'rows');\n",
    "display(edge);\n",
    "\n",
    "totalFace = uint32([elem(:,[2 3 4]); elem(:,[1 4 3]); ...\n",
    "                    elem(:,[1 2 4]); elem(:,[1 3 2])]);\n",
    "sortedTotalFace = sort(totalFace,2);                \n",
    "[face, i2, jf] = unique(sortedTotalFace,'rows');\n",
    "display(face);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In iFEM, `N,NE,NF,NT` represents the number of vertices, edges, faces and tetrahedrons, resprectively.\n",
    "\n",
    "    N = size(node,1); NT = size(elem,1); NF = size(face,1); NE = size(edge,1);\n",
    "\n",
    "In the assembling procedure, the matrix is always computed elementwise and then assemble to a big one. A pointer from the local index of a simplex to its global index is thus indispensible.\n",
    "\n",
    "**Elementwise pointers**\n",
    "\n",
    "- `elem2node = elem`\n",
    "- `elem2face(1:NT, 1:4)`\n",
    "- `elem2edge(1:NT, 1:6)`\n",
    "\n",
    "Such information is exactly stored in the output of `unique` function. For example, `elem2face(t,1) = 17` means the first face of t (spanned by `[2 3 4]`) is the 17-th element in the `face` matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem2edge =\n",
      "\n",
      "  6�6 uint32 matrix\n",
      "\n",
      "    1    2    6    8   10   12\n",
      "    3    2    6   11   13   12\n",
      "    4    5    6   15   16   18\n",
      "    4    7    6   17   16   19\n",
      "    1    5    6    9   10   18\n",
      "    3    7    6   14   13   19\n",
      "\n",
      "\n",
      "elem2face =\n",
      "\n",
      "  6�4 uint32 matrix\n",
      "\n",
      "   13    5    3    1\n",
      "   15    5    6    4\n",
      "   17   11    9    8\n",
      "   18   12    9   10\n",
      "   14   11    3    2\n",
      "   16   12    6    7\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N = size(node,1); NT = size(elem,1); NF = size(face,1); NE = size(edge,1);\n",
    "elem2edge = uint32(reshape(je,NT,6));\n",
    "elem2face = uint32(reshape(jf,NT,4));\n",
    "display(elem2edge);\n",
    "display(elem2face);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Face to edge Pointer**\n",
    "\n",
    "`face2edge(1:NF,1:3)` records the global indices of three edges of a\n",
    "face. This pointer depends on the ordering of vertices of faces and the\n",
    "indexing of local edges in a face. We list the following two important\n",
    "cases. Other combinations is possible but not attractive.\n",
    "\n",
    "- Ascend ordering.\n",
    "\n",
    "All local faces and local edges are ascend ordered.\n",
    "    \n",
    "    locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];\n",
    "    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];\n",
    "    edgeOfFace = [1 2; 1 3; 2 3];\n",
    "    locFace2edge = [4 5 6; 2 3 6; 1 3 5; 1 2 4];\n",
    "\n",
    "- Consistent ordering\n",
    "\n",
    "The local face is ordered such that the corresponding orientation is\n",
    "consistent with the induced orientation.\n",
    "\n",
    "    locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];\n",
    "    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];\n",
    "    edgeOfFace = [2 3; 3 1; 1 2];   \n",
    "    locFace2edge = [6 5 4; 6 2 3; 5 3 1; 4 1 2];\n",
    "\n",
    "The global one can be obtained from the composition of `elem2face` and\n",
    "`locFace2edge`. For example, for the asecnd ordering scheme,\n",
    "\n",
    "    face2edge(elem2face(:,1),:) = elem2edge(:,[4 5 6]);\n",
    "    face2edge(elem2face(:,2),:) = elem2edge(:,[2 3 6]);\n",
    "    face2edge(elem2face(:,3),:) = elem2edge(:,[1 3 5]);\n",
    "    face2edge(elem2face(:,4),:) = elem2edge(:,[1 2 4]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ordering of Vertices\n",
    "\n",
    "We discuss the ordering of vertices of simplexes. Again there are local\n",
    "ordering and global ordering. They may not be consistent and a sign array\n",
    "is used to record the inconsistency if any.\n",
    "\n",
    "The local ordering refers to the ordering of vertices in `locFace` or\n",
    "`locEdge`, i.e. the ordering of the local index of vertices. For elements\n",
    "associated to faces or edges, the local ordering could be used in the\n",
    "formulation of the local basis and thus the ordering does matter.\n",
    "\n",
    "The global ordering refers to the ordering of vertices in `face` or\n",
    "`edge`, i.e., the ordering of the global index of vertices. Note that\n",
    "that in either local or global ordering, permutation of vertices will\n",
    "represent the same simplex. To fix an ordering we need extra information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### elem\n",
    "The local ordering is always [1 2 3 4]. Any permutation of four vertices of a tetrahedon still represents the same tetrahedron. Such freedom provide a room to record more information like:\n",
    "\n",
    "* global ordering of vertices\n",
    "* an orientation of element\n",
    "* refinement rules (uniform refinement or bisection)\n",
    "\n",
    "For 2-D triangulations, three vertices of a triangle in 2-D is sorted counter-cloclwise and the first vertex is chosen as the newest vertex. Such ordering enables the efficient implementation of local refinement\n",
    "and coarsening in 2-D; see [Bisection in Two Dimensions](../afem/bisectdoc.html)\n",
    "and [Coarsening in Two Dimensions](../afem/coarsendoc.html).\n",
    "\n",
    "In 3-D, for the longest edge bisection, the newest vertex (with the highest generation) is stored as the last (4-th) vertex of a tetrahedron. For [3-D Red Refinement](uniformrefine3doc.html), the ordering determines the shape regularity of refined triangulation. Permuation of vertices in `elem` could deterioriate the angle condition after the refinement.\n",
    "\n",
    "We shall reserve the ordering of `elem` from the mesh refinement and\n",
    "coarsening since they are more subtle. We switch the ordering when\n",
    "generating data structure for finite element basis and assemble the\n",
    "matrix equation. Such sorting is hidden in the subroutines when a finite\n",
    "element basis requiring ordering is generated.\n",
    "\n",
    "Two types of ordering of `elem` is of particular importantance\n",
    "- Ascend ordering\n",
    "- Positive ordering\n",
    "\n",
    "In the ascend ordering, the vertices of `elem` is sorted such that \n",
    "\n",
    "    elem(t,1) < elem(t,2) < elem(t,3) < elem(t,4). \n",
    "    \n",
    "Such ordering will benefit the construction of local bases for high order basis or basis with orientation. This can be easily achieved by `elem = sort(elem,2)`. One has to rotate the boundary flag accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem =\n",
      "\n",
      "     1     2     3     7\n",
      "     1     3     4     7\n",
      "     1     5     6     7\n",
      "     1     5     7     8\n",
      "     1     2     6     7\n",
      "     1     4     7     8\n",
      "\n"
     ]
    }
   ],
   "source": [
    "bdFlag = setboundary3(node,elem,'Dirichlet');\n",
    "[elem,bdFlag] = sortelem3(elem,bdFlag);    \n",
    "display(elem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the positive ordering, the four vertices are ordered such that the\n",
    "signed volume, the mix product of vectors (v12,v13,v14), is positive.\n",
    "This is the default ordering used so far. `fixorder3` will switch the\n",
    "vertices for elements with negative volume. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem =\n",
      "\n",
      "     1     2     3     7\n",
      "     1     3     4     7\n",
      "     1     5     6     7\n",
      "     1     5     7     8\n",
      "     1     6     2     7\n",
      "     1     7     4     8\n",
      "\n"
     ]
    }
   ],
   "source": [
    "elem = fixorder3(node,elem)   % switchs the vertices for elements with negative volume."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### edge\n",
    "\n",
    "For 3-D triangulations, we chose the ascend ordering both locally and globally. Namely\n",
    "\n",
    "    locEdge(:,1) < locEdge(:,2); \n",
    "       edge(:,1) < edge(:,2);\n",
    "\n",
    "Recall that for `locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]`, it is ascend\n",
    "ordered. The `edge` produced by `unique` function is also ascend\n",
    "ordered. \n",
    "\n",
    "There might be inconsistency between the local and global ordering.\n",
    "That is `edge(elem2edge(t,1),1)` may not be smaller than `edge(elem2edge(t,1),2)`. It will be more clear from the discussion of the corresponding orientation.\n",
    "\n",
    "For 2-D triangulations, the global ordering is still ascend ordered. But locally it may not. For example, for the consisitent ordering `locEdge = [2 3; 3 1; 1 2]`, then `locEdge(2,1) > locEdge(2,2)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### face\n",
    "\n",
    "For 3-D triangulations, the `face` produced by `unique` function is already sorted in the second dimension such that the global ordering is ascended i.e. `face(:,1) < face(:,2) < face(:,3)`. The local ordering in `locFace`, however, is not always ascend ordered.\n",
    "\n",
    "    locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; % Ascend ordering\n",
    "    locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]; % Consistent ordering\n",
    "\n",
    "Again the local and the global ordering maynot be consisitent. That is\n",
    "\n",
    "    face(elem2face(t,:),1) < face(elem2face(t,:),2) < face(elem2face(t,:),3)\n",
    "    \n",
    "maynot be always true unless we use the ascend ordering in both `face` and `locFace`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orientation\n",
    "\n",
    "The orientation of a tetrahedron is either positive or negative. The\n",
    "orientation of a face is given by a normal vector and the orientation of\n",
    "an edge is determined by a tangential vector. \n",
    "\n",
    "The orientation of a d-simplex will induce an orientation of its d-1 boundary\n",
    "subcomplex and is called the *induced orientation*. For example, a positive\n",
    "orientated tetrahedron will induce the outwards normal orientation of its\n",
    "four faces and a positive orientated triangle will induce the counter\n",
    "clockwise orientation of its three edges.\n",
    "\n",
    "The ordering of vertices of a simplex will naturally introduce an\n",
    "orientation and will be called the *ordering orientation*. More specifically\n",
    "\n",
    "- the vector from `edge(:,1)` to `edge(:,2)` defines an orientation of edges.\n",
    "- the `cross(v12,v13)` defines an orientation of a face, where `vij` is\n",
    "the vector from vertex `face(:,i)` to `face(:,j)`.\n",
    "- the sign of the mix product `sign(v12, v13, v14)` defines an orientation for tetrahedrons.\n",
    "\n",
    "The orientation of a simplex in the simplicial complex should be uniquely\n",
    "determined which will be called the *global orientation*. It can be chosen as\n",
    "the global ordering orientation but not always.\n",
    "\n",
    "Inside one tetrahedron, the local ordering of local edges and local faces\n",
    "will introduce a corresponding orientation. The orientation of the\n",
    "tetrahedron will also induce an orientation for its four faces. These are\n",
    "called the *local orientation* which may not be consisitent with the global\n",
    "orientation. The local ordering orientation is used in the local basis and the induced orientation is used when computing the differential operator locally.\n",
    "\n",
    "In general, there will be an inconsistency of different types of orientation and apporipate data structure should be constructured to record such inconsistency.\n",
    "\n",
    "- a global orientation\n",
    "- the global ordering orientation\n",
    "- the local ordering orientation inside a tetrahedron\n",
    "- the local induced orientation inside a tetrahedron\n",
    "\n",
    "We now discuss the orientation of elem, face, and edge separately. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### elem\n",
    "\n",
    "The orientation of a tetraheron is either positive or negative.\n",
    "We chose the global ordering orientation, i.e., the sign of the signed\n",
    "volume computed from `elem`. \n",
    "\n",
    "    [Dlambda,volume,elemSign] = gradbasis3(node,elem);\n",
    "\n",
    "In the output of `gradbasis3`, `volume` is always positive and an\n",
    "additional array `elemSign` is used to record the sign of the signed\n",
    "volume.\n",
    "\n",
    "`Dlambda(t,:,k)` is the gradient of $\\lambda_k$ associated to vertex $k$. Therefore the outward normal direction of the kth face is `-Dlambda(t,:,k)` which is independent of the ordering and orientation of the element."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### face\n",
    "\n",
    "Again we use the global ordering orientation determined by `face`. The normal vector is given by `cross(v12,v13)`.\n",
    "\n",
    "The local ordering orientation is implicitly used when computing finite element basis in each element. For example, the RT0 basis on face `[i j k]` in `locFace` is defined as\n",
    "\n",
    "$$\\phi_{i,j,k} = 2(\\lambda_i \\nabla \\lambda_j \\times \\nabla \\lambda_k+ \\lambda_j \\nabla \\lambda_k \\times \\nabla \\lambda_i+\\lambda_k \\nabla \\lambda_i \\times \\nabla \\lambda_j).$$ \n",
    "\n",
    "Odd permutation of `[i j k]` will change the sign of the basis. The direction of $\\phi_{i,j,k}$ is the normal vector determined by `[i,j,k]` ordering. Note that this is defined locally, i.e., element by element. \n",
    "\n",
    "The global basis associated to a face, however, depends only on the global orientation of this face. We introduce `elem2faceSign(1:NT, 1:4)` to record the inconsistency of a local ordering orientation and a global orientation.\n",
    "\n",
    "For `locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]`, i.e. the induced orientation, the `elem2faceSign` can be obtained from subroutine `dof3face` follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem2faceSign =\n",
      "\n",
      "     1    -1     1    -1\n",
      "     1    -1     1    -1\n",
      "     1    -1     1    -1\n",
      "     1    -1     1    -1\n",
      "    -1    -1     1     1\n",
      "    -1    -1     1     1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "totalFace = [elem(:,[2 3 4]); elem(:,[1 4 3]); elem(:,[1 2 4]); elem(:,[1 3 2])];\n",
    "elem2faceSign = reshape(sum(sign(diff(totalFace(:,[1:3,1]),1,2)),2),NT,4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When both `elem` and `locFace` are ascend ordered, the orientation of the\n",
    "global ordering is consistent with that of the local ordering. Thus\n",
    "`elem2faceSign` is not needed for the ascending ordering in assembling the mass matrix.\n",
    "\n",
    "But for the asecond ordering system, an `elem2faceSign` will be used when assembling differential operators because the orientation for Stokes theorem is induced orientation. For example, when computing `div` operators on a positive orientated tetrahedron, the faces should be orientated by the outwards normal direction but the global faces may not be.\n",
    "\n",
    "If `elem` is positive ordered and `locFace` is consistently ordered, then this inconsistency is already recorded in `elem2faceSign`.\n",
    "\n",
    "For the ascend ordering of `elem` and `locFace`, we use $+1$ if the\n",
    "orientation of a face is the same with the induced outwords normal direction in a certain elem, and $-1$ otherwise. Then the inconsistency is given by `elem2faceSign = [1 -1 1 -1]` by comparing\n",
    "\n",
    "- The induced orientation: `locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]`;\n",
    "- The ascend orientation:    `locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]`.\n",
    "\n",
    "Here we use the *ascend orientation* to refer to the orientation given by the ascend ordering.\n",
    "\n",
    "In summary, \n",
    "\n",
    "- the induced orientation is favorable for computing $d\\phi$;\n",
    "- the asecond orientation is favorable for computing $(f, \\phi)$ or $(\\phi_i, \\phi_j)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### edge\n",
    "\n",
    "The orientation of edges is simpler than faces. Globally we always chose the global ascend ordering orientation. Namely the orientation of an edge is from the vertex with the smaller index to the larger one.\n",
    "\n",
    "Locally the local ascend ordering may not be consistent with the global one. See Lowest Order Edge Element. For the ascend ordering of `elem` and `locEdge`, the local and the global\n",
    "orientation will be consistent and no `elem2edgeSign` is needed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem2edgeSign =\n",
      "\n",
      "  6�6 int8 matrix\n",
      "\n",
      "    1    1    1    1    1    1\n",
      "    1    1    1    1    1    1\n",
      "    1    1    1    1    1    1\n",
      "    1    1    1    1    1    1\n",
      "    1    1    1   -1    1    1\n",
      "    1    1    1   -1    1    1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "totalEdge = uint32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...\n",
    "                    elem(:,[2 3]); elem(:,[2 4]); elem(:,[3 4])]);\n",
    "direction = ones(6*NT,1,'int8');\n",
    "idx = (totalEdge(:,1)>totalEdge(:,2));\n",
    "direction(idx) = -1;\n",
    "elem2edgeSign = reshape(direction,NT,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### face to edge\n",
    "\n",
    "For the ascend ordering `edgeofFace = [1 2; 1 3; 2 3]`, the local and the global ordering is consistent and so is the ordering orientation.\n",
    "\n",
    "Then it is not consisitent with the induced positive (counter clockwise) orientation of edges. When the edge direction is the same with the induced direction, we use sign $+1$, otherwise $-1$. Then `face2edgeSign = [+1 -1 +1]` records the inconsistency of the ascend orientation of the induced orientation.\n",
    "\n",
    "For the consistent ordering, `edgeofFace = [2 3; 3 1; 1 2]` which is consisent with the induced positive orientation but then may not be consistent with the global orientation of edges. We construct `face2edgeSign` to record such inconsistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "face2edgeSignp =\n",
      "\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "     1    -1     1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "totalEdge = [face(:,[2 3]); face(:,[3 1]); face(:,[1 2])];\n",
    "direction = ones(3*NF,1);\n",
    "idx = (totalEdge(:,1)>totalEdge(:,2));\n",
    "direction(idx) = -1;\n",
    "face2edgeSignp = reshape(direction,NF,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "We summarize the two popular ordering and orientation schemes below.\n",
    "\n",
    "### Ascend Odering and Orientation\n",
    "\n",
    "The asecond odering and orientation is more algebraic, determined by the indices of vertices. \n",
    "\n",
    "#### Ascend ordering\n",
    "The array `elem` is sorted such that\n",
    "\n",
    "    elem(i,1) < elem(i,2) < elem(i,3) < elem(i,4)\n",
    "\n",
    "The local face and local edges is also in the ascend ordering\n",
    "\n",
    "* `locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];`\n",
    "* `locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];`\n",
    "* `edgeofFace = [1 2; 1 3; 2 3];`\n",
    "\n",
    "Then due to the asecond ordering of `elem`, globally the `edge` and `face` also follow the ascend ordering, i.e.\n",
    "\n",
    "* `edge(e,1) < edge(e,2);` \n",
    "* `face(f,1) < face(f,2) < face(f,3).` \n",
    "\n",
    "One can easily see the benefit: the ordering of local edges and local faces is consistent with the global ones and so is their corresponding orientation.\n",
    "\n",
    "#### Orientation\n",
    "We chose the global ordering orientation for each elment. We chose the orientation corresponding to the ascend ordering for edges\n",
    "and faces. That is\n",
    "\n",
    "* `elem: sign(v12,v13,v14)`\n",
    "* `face:`  the normal vector is given by `cross(v12,v13)`\n",
    "* `edge:`  from the node with the smaller global index to the bigger one\n",
    "\n",
    "For faces and edges, the orientation of the ascend ordering and the induced orientation is not consistent. The inconsistency is recorded by\n",
    "\n",
    "* `elem2faceSign = [1 -1 1 -1];`\n",
    "* `face2edgeSign = [1 -1 1];`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Positive Ordering and Orientation\n",
    "\n",
    "The positive orientation and ordering is more geometrically consistent in the sense that the orientation of an element is locally consistent with the orientation of the local boundary faces. But it introduces inconsistency with the global orientation of a simplex. \n",
    "\n",
    "#### Positive and consistent ordering\n",
    "\n",
    "The vertices of `elem` is sorted such that the signed volume is always positive, i.e. the four vertices follows the right hand rule.\n",
    "\n",
    "The four faces of a tetrahedron are ordered consistently as \n",
    "\n",
    "    locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];\n",
    "\n",
    "The six edges of a tetrahedron still ascend ordered\n",
    "\n",
    "    locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];\n",
    "\n",
    "Three edges of a face is ordered consistently\n",
    "\n",
    "    edgeofFace = [2 3; 3 1; 1 2];\n",
    "\n",
    "#### Orientation\n",
    "\n",
    "The ascend ordering orientation is used for the global orientation of `edge` and `face` arrays. The inconsistency of the local and the global orientation is recorded in\n",
    "`elem2faceSign` and `elem2edgeSign`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Example \n",
    "We show two tetrahedron with the ascend ordering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem =\n",
      "\n",
      "     1     4     5     8\n",
      "     1     4     5     7\n",
      "\n",
      "\n",
      "edge =\n",
      "\n",
      "  9�2 uint32 matrix\n",
      "\n",
      "   1   4\n",
      "   1   5\n",
      "   1   7\n",
      "   1   8\n",
      "   4   5\n",
      "   4   7\n",
      "   4   8\n",
      "   5   7\n",
      "   5   8\n",
      "\n",
      "\n",
      "face =\n",
      "\n",
      "  7�3 uint32 matrix\n",
      "\n",
      "   1   4   5\n",
      "   1   4   7\n",
      "   1   4   8\n",
      "   1   5   7\n",
      "   1   5   8\n",
      "   4   5   7\n",
      "   4   5   8\n",
      "\n",
      "\n",
      "elem2edge =\n",
      "\n",
      "  2�6 uint32 matrix\n",
      "\n",
      "   1   2   4   5   7   9\n",
      "   1   2   3   5   6   8\n",
      "\n",
      "\n",
      "elem2face =\n",
      "\n",
      "  2�4 uint32 matrix\n",
      "\n",
      "   7   5   3   1\n",
      "   6   4   2   1\n",
      "\n",
      "\n",
      "v =\n",
      "\n",
      "    0.6667\n",
      "    0.6667\n",
      "\n",
      "\n",
      "elemSign =\n",
      "\n",
      "    -1\n",
      "     1\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4wEFCiILa00AcAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNS1KYW4tMjAxOSAwMjozNDoxMYk69nwAABe2\nSURBVHic7d1vjF1lnQfw506HLhUnMpJOdWsRS7e4BsKfJgqyAYeFaOoWibLiv2AnxKDBVYbF0LRk\nVxSx3bo7YOLuClYrY1oETDBdaxe6dt/AUuBFQQELTSPyTzs2g5Q/pUw7++KU6eXOnTv33rn3nuec\n8/m8ur233HnelO98f89zzimNj48HAEhbV9oLAIAQBBIAkRBIAERBIAEQBYEEQBQEEgBREEgAREEg\nARAFgQRAFAQSAFEQSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRA\nFAQSAFEQSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQ\nSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQSFCXAwcO\nPP/882mvAvKsO+0FQDZs3bp1YGDgvPPOe8tb3jJnzpzyj6666qqFCxemtTDIDYEEddm7d++ePXtu\nu+22ivdPOeWUG2+8MZUlQc4IJKjL9u3bf/KTn7zrXe/q6joy6L766quHh4e7u/07ghbwDwnq8uqr\nr372s58tf2fdunUDAwOLFy9Oa0mQMwIJ6rJq1aryP+7atWvr1q0bN25Maz2QP6Xx8fG01wDZ85GP\nfGR4eHju3LlpLwTyw7FvaNjNN9/8gQ98QBpBa2lI0JhDhw69+93vvueee9773vemvRbIFXtIMKVS\nqTTxeuJXt40bN46MjEgjaDmBBFWUR1H5O+Pj4z/96U9dBgvtYA8JGnDo0KEtW7YsWrQo7YVADmlI\nUGlyPZowa9asEMLs2bM7uBwoCg0JGnbUUUelvQTIIYEEDXOvIGgH/66gMRs2bDjnnHPSXgXkkOuQ\noIqKbaS+0/r27NiTvPZPBtrEyA7q0ndaXwhh89Obaxx5AGbCyA6qSGpQefas37Q+eTGRSaoStJaR\nHUyjVCptfnrz5PeXLljqnw+0kJEd1FIqlZJh3WRJVTLBg1YxsoPmJc2pVDJpgBbQkGCmnHSAlhBI\n0ALGdzBzRnbQGsZ3MEMaErSSqgRNE0gwpakOfNe2+enNdpWgCQIJ2kJVgkbZQ4J2sasEDdGQoL1U\nJaiThgRtpypBPTQk6BBVCWrTkKC65o7Y1aYqQQ0aEnSaqgRVaUiQAlUJJtOQIDWqEpTTkCBNqhJM\n0JAgfW41BEEgQVXtOGJXm/EdGNlBLIzvKDgNCeKiKlFYGhJER1WimDQkiJSqRNFoSBAvVYlC0ZCg\nUueP2NWmKlEQGhJkgKpEEWhIkBmqEvmmIUGWqErkmIYE2aMqkUsaEmSSqkT+aEjwJrEdsavNXVnJ\nE4EE2WZ8R24Y2UHmGd+RDxoS5ISqRNZpSJAfqhKZpiFB3qhKZJSGBEdk64hdDaoSWaQhQW6pSmSL\nhgR5piqRIRoS5J+qRCZoSFAIqhLx05CgQFQlYqYhQbGoSkRLQ4LDcnPmux7uykqEBBIUlPEdsTGy\ng+IyviMqGhIUnapEJDQkQFUiChoScJiqRLo0JAihYEfsalCVSJGGBFRSlUiFhgRUoSrReRoSMCVV\niU7SkIBaVCU6RkMCpqcq0QEaEjhiVxdViXbTkIAGuCsr7SOQgMZsfnrzvNPniSVazsgOaNj6TeuD\n8R2tpiEBTXLSgdbSkIDmOelAC2lIFJ0jdjOnKtESGhLQAqoSM6chAS2jKjETGhLQSqoSTdOQgNZT\nlWiChgS0hapEozQkCs0Ru3ZTlaifhgS0l6pEnTQkoBNUJaalIQEdoipRm4YEdJQHWDAVgQTNeP75\n539wyw/SXkVWGd9RlZEdxTWTI3bDw8M9PT2tXU+hGN8xmYYEDbv77rvHxsbSXkUeqEqU05BgSoOD\ngxMvhoaGktfPPfdcV1dX39y+18deT29p+aEqMUFDgioGBwcn0qjinV9t+9X555+f0rpyS1UiaEjQ\nkLvvvvuCv70g7VXkk6qEhgSVKrpRuV/+8pfz3jGvk4spGlWpyDQkCspd7KKlKhWWQILpje4aTXsJ\nhTNx/axYKg6BBNPrXdQ78Xr5suUTr7dv3l77P1y/aX17VlQITVelXbt2feELX/j6179+7rnntmdp\ntIVAgkpDQ0MV20g983v2PbuvZ37PyaecnLzz+GOPz5kz54T3nHDUUUd98KwPVnzDvDmH95nK02ta\n0quq+qvSE088cd111+3Zs+ell166//77X3/dufyMEUjQgHPPOfwb964ndx177LETf6zwx1f/mLxY\nc/ua5MU1n7xm8l+b+DQJsMHBwfpng4VKrzqr0qJFi3784x93d3cPDw/ff//9nVodLSOQoIrkMtgk\nIZJ53YobVqxeuXriLxw6dOi1116b9nvKc2jbtjd91N9/5NMkmVbcsKKetRU2vZKqVJFJ5efxxsfH\nu7qcHM4w51goqDpP2S1ftnyix4QQVq9cffqHTv/TyJ9eeOGFEEJvb+/cvrkXnF/9yqRrPnlNRQhN\npb8/lP+UFmoivULcAbZ0wdIQwvj4+FRHw8fHx4eHhy+99NJ77rnHJczZoiFRRE2f+U56Uj1Vpv40\nCiFs2xb6+69pRyYlw8OGuleIe+trYnzXyR9KZwgkaL2G0ijRvkyq3+StrxrSTa+TLj5p5507q35U\nKpVuvfXWlvwUOkwgwZQq5nWJaUtSE2mUiCGT6tdQeoU3AqxV6TW6a7TvtL7eRb1TxRJZJJCgYTUy\nqek0SmQrkxqSBFhD6TXt1lcSS3t27GnJCkmdQIJm9Myv8nS+GaZRIseZVL9pt77KTzwKpNxwRBKa\nccUVV5T/P7FRBw+GX/wirF0bvvrV8K3rw6ZN4aAH/jVixQ0reub39Mzv2XnnzpMuPqn8IyeHs0tD\nonCaO/A9Wf0n7iq8MBr+6Z/Dr3995J2t/xN+9rNw003hbW9r9MsK7aGbHqpIIzJNQ4LmlV8tW/Ve\nDFV9/+bDadTdHfr6Dr/51FNh7b8c+Tv1f1sxDSweeOimhyp+sRgfH0/q0b59+0IIHjOfORoStEw9\nG0gjI+Hu/w4hhLceE7737+H448P/3RdWrgohhO0PhAMHwuzZyTZSe5eaaQOLByaiKLmnxtI7lyZR\ndNFFFz3zzDM7duwIIVx44YVnn332woUL161bl+JqqZ9AgioGBwfXb1o/cbK5hkYHd79/KhwaDyGE\n/vPC8ceHEMJZHwzz/zI8+1wYGwvPPRdOOKHpVRdCeRpNdtddd3VyMbSWQIIqRneN1pNGiYrb3NX2\n+lhYtCiEEE499fA7+/eHP/85hBC6SuGd72hwoUUysHggvHGnhnJLFyx1kCEfBBK0wIobVtS563Pm\nmeHMM4/88cBr4VvfCi+9HEIIp58R/uLo9qwv+2oXI/JBIFEsUT25fNeT4durw+7dIYRw7LHhH69K\ne0FRmqoYkT8CCSotX7a8zg2kCv39dZ1rSNx2W1i3LiQHwRYtCt/8RnjHO498DwnFqFAc+4Yqmkij\n+u+tcHAsfOO68P3vh7GxMKsrfP7z4T//40gaNfptOSaNikZDgk77138L2/43hBB6esINN4STT055\nPREypismgQQd9fjjYcuWw69PXBjuuzfcd++RTz92UZg3L5V1RUQxKiyBBG/S9AZSCGHN7Wv6+6e5\nv+r27WHiiPKOh8OOh9/06RlLwqc+Veh5nTQqMntIFEidR+yaS6PEmtvX1D6S8MwztT792teKm0YD\niweaSCMXIeWJhgQtVrsnXXttuPba6h/19xc6jRQjNCRovWl70mTSKO1VkD6BBEdM+8iJ+jWUSYVN\no+bGdOSVkR20SzK7m/hjxRCvPK4Km0aiiHICCdpoze1rknuBX/PJayYXpmLmUEIaMZlAoijSuotd\nci/wImdPBRe9MhV7SHBYCzeQKjT0fIp8S4qRNKIqgQSdIJOCMR3TMbID2q5NYzpXxeaMhgQhvPHM\n8rb+iMKWJGM66iSQIIQGn1netKJlksuMaIhAohDieVBscTJJMaJRAgloPcWIJggkOPzIiY79uHyX\nJGM6miaQIISZPXKiCXnNJGM6ZkIgQTryl0mKETPkOiRgplK5G5CLkPJHQyL/ah+x6/AGUrl8lCRj\nOlpFIEGnN5DKZT2TjOloISM7oBlu2k3LaUiQsiyWJGM62kEgUWjte+REQ7KVScZ0tIlAgihkIpNc\n9Epb2UMC6iKKaDcNiZyL57aq04q5JEkjOqDkyjLyrUYgRbKBVGH1ytUrbliR9iqOiPM0natic0lD\ngrhE1ZOcpqOT7CEBVcRZjMg3DYmC6sAzy5uWeklSjEiFQKKgOvPM8qalmEnOL5AWgUSeZeiI3WSd\nzySXGZEue0hACIoREdCQKKIUHznRkI6VJGlEDAQSBRXzBlK5dmdSFsd0LkLKKyM7KK7MRRH5piFB\n7NpUkqQRsRFI5NZUR+yysoFUrrWZlMUxHUUgkCiirGwglWtVJrnolWgJJCgQxYiYCSTIjJmUJGM6\n4ieQKJY4HzlRv+YyyZiOTHDsGzImyaSJZyZVzafyJyopRmSFQCKfMn0Xu3qU59DQ0Js+Ghw8/OnO\nO3eG3D1CwlWxOSaQIJMqQmjyR0sX7Hx0tK+/v2Mrgpmyh0SBZH0DKbF65eoaaZRI0iiEsG1bWL5s\neQdWBTMnkCBLpk2jpQt2TqRRQiaRFQIJMqOeNHp0tK88jRIyiUwQSBRFzM8sr0ftNJpcjCrIJOIn\nkMihqkfsIn9meW3TplHVYlRBJhE5p+wg22oUow0bwquvhNNOD0uWdHhR0AyBBFm1dMHOEMJUafTb\nx8Mtt4QQwmsH8hNILkLKNyM7CiGLj5yYUPVeDDXGdAcPhvvuC6uurf5tpnZES0OiKLK7gRQmXQZb\nY0y3+tth27Zw4PXq37NtW3CpLNESSJAltcd0IYQXX5wyjSByAom8yfFd7Gof7E5c/Pfh3A+FEMKj\nvwmb/qsDi4KWEUjkX7KBlOmRXagvjUIIZ5xx+EVXl0AiYxxqoBAynUY775zmolfIBw0JojaweKDv\ntJZFkRMNxExDgnglz9Zr7YH17B5/J/c0JHIuo4+cGFg8EHL3bL0ZclVs7gkkciUfR+w8dJxiMrKD\nuFRNo/Wb1s98+6e/37yOqGlIEIvaY7r1m9b39y/ftq3eb+t682+b0oj4lcxkyZOKkV2GNpDqHNMt\nX9ZAJk3IRxrZQ8o9DQlS1tD5hUZ7UnDUm+ywhwRpSopRQ0cYGtpPSrrR+k3r3eSb+BnZkSvlI7vB\nwcGhoaGY79Ewk9N05QFTUZjK46p8UpfpZ3AEI7sCEEjkR4Y2kFp4mVHV6jNV8GQ3k6RREdhDgk5r\n7WVGGQ0YmMweEnRUuhe92kwiZgKJfIpwNjWweCCGWzDIJKIlkMitqI4zNHGarn1kEnESSNB2MRQj\niJ9AIifivK1qJGO6yZQkIiSQyKFINpCiGtNNJpOIjUAin1LfQIqzGFWQSUTFdUjQYp6t13Kuii0I\nDQlaKfIx3WRKEvEQSORNincMysSYbjKZRCSM7MiD1I/YZX1M17uoN+0lgIYEM5a5Md1kQ0NDShKp\nE0jQvGgvM2qCwR2pE0jkSic3kHJQjCrIJNIlkKAZuSlGEA+BBI3J05husghLkouQikMgkXkTR+wG\nBwfbfceg/I3pJoswkygIgUR+jO4abesdg3JcjCrIJFLhOiSYXtYvM4JM0JBgGkUY002mJNF5Aomc\naNMjJ4ozpptMJtFhAon8aO0GUr5P09VJJtFJ9pDItjbdxU4UQedpSFBJGpVTkugYgUQetGoDyZiu\nqhQzyVWxhWJkR07MfANJFNXg+RR0gIYEIUij6Xg+BR0gkCg6Y7o62Uyi3QQSGZYcsZvJIyeKedFr\n02QSbWUPiYJyNyCIjYZEESlGTVOSaB+BRLY18cgJO0YzJJNoE4FEtjX0yAnnF1qlM5nkIqSisYdE\nUYgiiJyGRFY1dBc7adRyBne0XEkjJqNKpdJJF580NDRUe2TnNF1btempHwkju6LRkMiwaTeQnKZr\nNz2JFhJI5JYxHWSLkR1ZVSqVuo/uPu/T55146onds7v3v7T/N/f+5j2nvOeCSy8wpuuwNg3ujOyK\nRiCRVaVSac5b57z60qvJH7tnd1/0Dxd99PKPKkapaEcmCaSiceybDDv+r48/95Jz97+8/+hjjj75\nb07ueXuPNEqL51MwcxoSGVMqlSZeL7lgyZe/9+XktTFd6lpbktSjAhJIZEZ5FJX70RM/Uowi0cJM\nEkgFJJDIjKkCKYTwpW9+6Xe//d3HLvvY8X91fCeXxGStyiSBVECOfZMNNdIohLBs+bIPf/rDV/7d\nlY89+FjHlgS0lkAiJxafuvjUs0/9zpXfSXshRedSWZomkMiPt/e9/Q+//8Mj9z2S9kKKTibRHMe+\nyYOlC5ZOvN77x70proREkkntu80duSSQyIbx8fGptpEmtr4vu+yyH/7wh2u/snbtV9aW/wUH8CAT\nBBL5MTIyMmvWrGeffXbevHnl709OMhHVATMpSY7YFZNAIjOS/0OVp8vPf/7zCy+8MHk9MjKyZcuW\nVatWVaRRKKtQE0RUZxjc0RDXIZFhV155ZW9v71lnnXXgwIE1a9acffbZq1evbu6rRFT7NJFJGlIx\naUhk2I033vjII488/PDD3d3dGzZsWLBgQdNfVU+LenT00adeeqrpHwHUpiFBvaqeqlCk6tFoSdKQ\nikkgwYyY9dWpoUwSSMVkZAcz4sREnTyfgmlpSNB2IipRf0nSkIpJIEEKChtRdWaSQComgQRRKM6h\nvmkzSRoVlj0kiEI9e1GhMEWKYtKQIEvyMeurXZI0pMISSJBtGY2oGpkkkApLIEHeZCWipsokgVRY\n9pAgb1waRUZpSFBEkRzqq1qSNKTCEkhACOkd6pucSQKpsAQSMKXOzPrKM0kaFZk9JGBKndyOWrpg\n6cSP2Ldv3/XXXz86OnrgwIGXX375E5/4xCWXXNKSn0LMNCRgRmYeUcuXLU9e7Nmxp/z95P9OY2Nj\nn/nMZ97//vdfffXVza+SLBBIQIs1GlG1AymEMDIycuKJJ7744ostWyJREkhA29U41Dc4OJi8M7pr\ntCKQwhuZdOjQodmzZ4+NjbV/paTJHhLQdnXeqa/vtL7JmRRC2Lhxoz2kIhBIQAomIqo8mSan0dat\nW0MId9xxx6233tqxtZEWIzsgTVWrUrmFCxc+8MADxx13XGfWQ4q60l4AUGg1fideu3Ztb2/v7t27\n3/e+9917772dXBWp0JCAlE1VksbHx0dGRq644oo77rijt7f3ySef1JPyTUMCUjY+Pl7xm/GmTZuS\nd+bOnXv77bdv2LBhdHR0eHg4pQXSIRoSEJFXXnnlmGOOOXjwYFfXm35d/vjHP37cccfdcsstaS2M\nDtCQgIgcffTRs2bN2rt3b8X7vb29S5YsSWVJdIxAAiLS1dU1MDCwcuXK8jf37t372GOPfe5zn0tr\nVXSGkR0Ql/37919++eV79+794he/OHv27N27d991113f/e53Fy9enPbSaC+BBMRo9+7dDz744P79\n++fPn3/++eenvRw6QSABEAV7SABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEg\nARAFgQRAFAQSAFEQSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRA\nFAQSAFEQSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQ\nSABEQSABEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQSABEQSAB\nEAWBBEAUBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQSABEQSABEAWBBEAU\nBBIAURBIAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQSABEQSABEAWBBEAUBBIAURBI\nAERBIAEQBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFEQSABEQSABEAWBBEAUBBIAURBIAERBIAEQ\nBYEEQBQEEgBREEgAREEgARAFgQRAFAQSAFH4fwLPP2uQ2CHGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "% A mesh with two tetrahedron with the ascend ordering\n",
    "elem = [1 4 5 8; 1 4 5 7];\n",
    "node = [1,0,0; 1,1,1; 1,-1,-1; 0,1,0; -2,-1,0; 1,1,-1; 0,1,1; 0,-1,-1];\n",
    "NT = size(elem,1);\n",
    "showmesh3(node,elem,[],'FaceAlpha',0.25);\n",
    "findelem3(node,elem);\n",
    "findnode3(node,elem(:));\n",
    "display(elem);\n",
    "% generate edge array\n",
    "totalEdge = uint32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...\n",
    "                    elem(:,[2 3]); elem(:,[2 4]); elem(:,[3 4])]);\n",
    "sortedTotalEdge = sort(totalEdge,2);\n",
    "[edge, ~, je] = unique(sortedTotalEdge,'rows');\n",
    "display(edge);\n",
    "% generate face array\n",
    "totalFace = uint32([elem(:,[2 3 4]); elem(:,[1 4 3]); ...\n",
    "                    elem(:,[1 2 4]); elem(:,[1 3 2])]);\n",
    "sortedTotalFace = sort(totalFace,2);                \n",
    "[face, i2, jf] = unique(sortedTotalFace,'rows');\n",
    "display(face);\n",
    "% generate pointers of indices\n",
    "elem2edge = uint32(reshape(je,NT,6))\n",
    "elem2face = uint32(reshape(jf,NT,4))\n",
    "% find orientation of elem\n",
    "[v,elemSign] = simplexvolume(node,elem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are using the ascend ordering, the inconsistency with the induced orientation is \n",
    "* `elem2faceSign = [1 -1 1 -1];`\n",
    "* `face2edgeSign = [1 -1 1];`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary Faces and Boundary Conditions\n",
    "\n",
    "We use `bdFlag` to record the boundary condition; see [Data Structure: Boundary Conditions](bddoc.html) for details. For short, `bdFlag` has the same size with `elem`, and records the boundary type of each local faces. If we change the ordering of `elem`, the corresponding local faces are changed. Threfore when we sort the `elem`, we should sort the `bdFlag` respectively. We use `sortelem3` to sort `elem` and `bdFlag` at the same time. Note that `sort(elem,2)` sorts the `elem` only, and leave\n",
    "`bdFlag` unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "elem =\n",
      "\n",
      "     1     2     3     7\n",
      "     1     4     3     7\n",
      "     1     5     6     7\n",
      "     1     5     8     7\n",
      "     1     2     6     7\n",
      "     1     4     8     7\n",
      "\n",
      "\n",
      "bdFlag =\n",
      "\n",
      "  6�4 uint8 matrix\n",
      "\n",
      "   1   0   0   2\n",
      "   2   0   0   2\n",
      "   2   0   0   2\n",
      "   2   0   0   2\n",
      "   1   0   0   2\n",
      "   2   0   0   2\n",
      "\n",
      "change to ascend ordering\n",
      "\n",
      "elem =\n",
      "\n",
      "     1     2     3     7\n",
      "     1     3     4     7\n",
      "     1     5     6     7\n",
      "     1     5     7     8\n",
      "     1     2     6     7\n",
      "     1     4     7     8\n",
      "\n",
      "\n",
      "bdFlag =\n",
      "\n",
      "  6�4 uint8 matrix\n",
      "\n",
      "   1   0   0   2\n",
      "   2   0   0   2\n",
      "   2   0   0   2\n",
      "   2   0   2   0\n",
      "   1   0   0   2\n",
      "   2   0   2   0\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4wEFCiU2fGTapgAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNS1KYW4tMjAxOSAwMjozNzo1NI+JaBkAACAA\nSURBVHic7d1/jJXVve/xNfwYAcHKL5kWWypSy9D6AyeRTCW59xCqBiEaSY6XHz2t5zQ0kaREwzXn\nHrpjuVMNsX/gGKZJNf4jtKO19Ec6NC1M9bbVcKCdSOEc7Kn7AgrUkVHBO8hQhmHfP5768Lj3nr2f\n/fxY67vWer/iH8OeYc9SZ/Znf9Zaz7OaSqWSAgDAtDGmBwAAgFIEEgBACAIJACACgQQAEIFAAgCI\nQCABAEQgkAAAIhBIAAARCCQAgAgEEgBABAIJACACgQQAEIFAAgCIQCABAEQgkAAAIhBIAAARCCRA\nugsXLrz99tumRwHkbpzpAQCoo7e394EHHliyZMmkSZMmTpwY/dTDDz88d+5cUwMDskUgAdK99957\np06dev7558sev/HGG5988kkjQwLyQCAB0u3bt2/Hjh3XXnvtmDGX59g3bty4ffv2ceP4FYY7+GkG\npBsaGlqzZk30kWefffaBBx644YYbTA0JyAOBBEi3adOm6B+LxWJvb293d7ep8QA5aSqVSqbHAKAB\nd9111/bt22fOnGl6IEDG2PYN2OTpp59etGgRaQQn0ZAAa1y6dGnOnDl79uyZP3++6bEA2WMNCRCq\nqakp/Dh449jd3T0wMEAawVUEEiBONIqij6xYsYLLYOEw1pAAa/zqV7+aN2+e6VEAeSGQAFkq61Fo\neHi4ublZ52AAnQgkwCbjx483PQQgLwQSYBPuFQSHEUiAUHd87w6l1Lzl89r/rT18cMuWLeZGBOSL\nd1uACOHSUaGvsG/fPqXU4IlBpdTcZXMHTwze8b07dj+4Wyk1e/Zsg4MEckUgAYYFUdR1uKt/qD98\ncPDE4JRrp0S/7I7v3XHPf7+nqYmL2eEsAgkwI1qJlFLRNFq0aFHvid7wj1OunTJ4YnDpyqX9Q/2F\nvkLwF4kluIdAAnSrWomienf23rbits72znnL5ymlWhe27j+xP/xsEGBUJbiHQAI0qVGJooI0Knvw\nthW39e7sXbpyafgIVQnuIZCAfJXlUIaoSnAMgQTkpe7UXKWwHnW2d657Zd1LW14KP1VZkgJUJTiD\nQAIyFnNqrlLVybo4wqqkiCXYjEACMhNEQlZTc4d+eWjJvy45dvBY+MhoJSlQ6Cu0TGxhBg/2IpCA\ntDJZJaqsRwMHB9Syxp6EfeGwGoEEJJThboWyNKpcQArVLknheKhKsBGBBDQswW4FzahKsBGBBMSV\neLdCbY3uZYhTkgJUJdiFQALqyO9CIpViZ11MVCVYhEACRmVkai5YQKr9NfFLUoCqBCsQSEC5nKbm\nKuVdj6KoSpCPQAIu01mJUqZRoyUpwN2GIBmBBOirRKG6afTSlpeW/OuSnL47VQkyEUjwV667FZKJ\ns4AUSlaSAlQlCEQgwUdmLyTSuXRUG1UJohBI8Ij+qblK2aZRmpIUoCpBDgIJXpB/bwWzqEqQgECC\nyyRUoqja9aihBaRQ+pIU4AwLGEcgwUECdysoSUtHNXAJLQwikOAUpubS4xJamML7ILhAZiWKqluP\novN1h3556MZlNwYfHzt4rHVha93n3/+L/eln7cq0TGxZv2A9LxHQhoYEu1lRiRqdrEtwNF8eqErQ\njECClaTtVjAuq60NlVhVgjYEEmwif2qukhV7GWqjKkEPAgl2sGJqrlLMNEq24btMfiUpQFVC3ggk\niMbUnChUJeSKQIJQllaiKP2TdXmXpAB3G0JOCCTI4kwlcmDpqDaqEjJHIEEEG3crZCWTBaSQnpIU\noCohWwQSDHNgaq5SmnqU69F8eaAqISsEEsxwZmqukvHJOp0lKcCNWZEJAgm6BS9bHk7NOY994UiJ\nQIImnqwSNVqPsl1ACukvSQH2hSMNAgn58iSHAsYn64SgKiEZAgl5cXK3gkVMlaQAVQkJEEjImMO7\nFWqjHlWiKqEhBBIy4/NuhWRplNMCUshsSQpQlRAfgYS0vFolytuhXx6y6yKkmKhKiINAQkLkUCjD\nybrMj+aTUJICVCXURSChYexWiEqcRnnP18lEVUINBBLi8na3gu3klKQAVQmjIZBQB1NzNbCzLjFu\nzIpKBBJGxdRcbRalkbSSFKIqIYpAQjmm5uJImUZ+LiBVRVVCiEDCZVQih4ktSQGqEhSBBEUlapxF\nk3UW4QwLEEiinT179jvf+c7p06cvXLjw4Ycfrly58v7778/qydmtkEx+aWTd0Xx5YF+4zwgkuS5d\nuvS1r32ts7Nz9uzZSqmLFy+uXr36+PHjGzduTPnMTM2ZZWoBSfisXYh94d4ikOTauXNnW1tbkEZK\nqXHjxnV1dV1//fWJA4mpufSYrNOGquQhAkmu3/72t6dOnYo+Mn369HPnziV4KipRJhxII1tKUoCq\n5BsCSa729va1a9d+85vffOKJJyZMmKCU6u7ubmgNiUoEB1CV/MH/Y7kuXLjwxS9+8Y033vjsZz/7\nzDPPKKW2bdv23HPPXXXVVbX/IrsV8pBVPaqxgFR1U8Oxg8daF7am/75R+3+x35aSFNXR1qGoSk6j\nIcnV3Nz8+9//fu3atb29vV/+8pfnzp27f//+2mnE1FxOHJiscwBVyXljTA8AtbzzzjtTpkz57ne/\nO3Xq1CNHjixYsODVV1+t/LKmjxT6CoW+AmmE2oKVJNOjSCJcVQqnAeASGpI40d+0m2+++Te/+c30\n6dO/+tWvrl+//sUXX1yxYsUbb7wxffr06BdTiXKlpx65ejRfHqhKrqIhCVL5vu9Pf/rTjBkzlFIz\nZ8780Y9+9MMf/vD06dPbt2+nEmmTbRrVWEAaODiQ1XeJw96SFKAqOYmGZJPVq1crpR566CF2KwCK\nG7M6h4YkRY03emEf6jrcdcs9tyz7X5mecY3Rub2XwfaSFKIqOYOGZIfgneDR/qMDRwbu3Hin6eF4\nIfM04siJ/FCV3EAg2eHIvx85feL0f/2f/7pn8z3Nk5pNDweOsOvGDXVxWwfbEUjSzVky529n/jb4\n7uC0T09bvW216eH4wu3JOodxhoXVCCQpSqVS5ST4vOXzij1FtjBo5lUaOVaSAuwLtxSbGqSbt3xe\nR1tHcNMU2IsFJM3YF24jGpIgwbu56O/P3+vRZqWU6mjroCppYKQemT2az8mSFKAq2YWGJE6pVAp+\neboOd0UfL/QVqEp582qyzhNUJYsQSKJ1He6KJlBwXwYyCZlz5pqk0RT6Cl2Hu8gk4Qgk+1CVcpJf\nPWIBSQKqknwEknRlJSlAVcock3XOl6QAVUkyAsliVCUgAaqSWASSBaqWpABVKRPUo4AnJSkQ/O6Q\nSaIQSC6gKqWRdxqxgCQZVUkUrkOyQ9fhrvUL1te4Din4VJBJXK4Un4RuJOpoPoevSRoNN2aVg4bk\nFGbwbKT5aD5URVWSgECyRo2VpDLM4MUkoR4J5NVKUlS4qkQsmcKUnZvCGTym70ajJ41YQLIOdxsy\niIZkk/glKUBVAhJgX7gpBJLjWFWqism62rydtYviElr9CCTLNFqSAlSlKG1pxHyd7ahKmhFIvqAq\noSGUpBBVSRsCyT7JSlKAqsRkHRKgKunBLjvv+LwBT2AamT2arzYPL5Idjc/v4XQikKxU98YNdYXT\ndx7Gkh4sILkh+mtCLOWNQPKXb1VJYD2Sz9uSFGZP+NvR0dYRvBE0Nyj3EUi2Sl+SAp5UJdIIMVX9\ndQjSyNCIPEIggRuzohZPSlJlJYJ+BJLFug53bX1k66rNqzJ5Nodn8PTXIxaQbBEnh6hH2hBIdiv2\nFNXmLJ/QvRk8JuvSc7Ikxfw5J410IpBQzuGqBDA1JxmBZLestjZUcqMqCa9Hoo7mq82BkpTg55l6\npBmBhFHZXpVMpVH8BaSBgwNqWd7D8V3iSkQa6UcgWS+/khRwoyohPbtKElNzNiKQUJ+NVUn4ZB3y\nk8n7J+qREQSSC/IuSQGLqhJplBPJJSnDSkQamUIgoQE2ViXNuAJJM6bmXEIgOUJPSQoIr0rUo1zJ\nKUk5/RBSjwwikJCE2KpEGjkv10pEGplFILlDZ0kKCK9K+nkyX2eqJPHD5jwCCamIqkp21SPJR/OJ\nom2ViHpkHIHkFP0lKSChKtmVRrbTUJI071YgjSQgkJANzrBAVvgp8haB5BpTJSlgagZPQj3yZAEp\nlHlJMriBm3okBIGE7GmewZOQRkjDbCUijeQgkBxktiQFtFUl0sheXNOKMgQSciRhswPyk2zWTlQO\nUY9EIZDcJKEkBXKtSnLqkW8LSMlIe3dCGklDIEGHPKqSnDRKwKKj+WqLU5JEVSJIRiA5S05JCoi6\nhNY4T47mk1aJoqhHAhFI0CqrqmR1PXJMZUmSX4lII5kIJJdJK0mB9FVJWhqxgBSQn0MQjkCCGWzA\nc8nex/fufXyvLf8rqUdiEUiOk1mSAsmqkrR65LPO9s7ggw17N+z/xX6zg4mJNJKMQIJhDVUlgWnk\n53xdEEUb9m4IH5FzcF8N3Y92k0aSEUjuk1ySAmzAs0W0EpkdSTLFnqJ6wvQgMDoCCVLUrUoC65En\nYuaQ8JLEZJ18BJIXug53bX1k66rNq0wPpI4aZ1i4lEYWHc1XOTVnKdLICgSSL4o9RbXZ9CDisWgG\nz9UFpMRTc8JLEoQjkCBUdAbPpXoknDOVKIp6ZAsCyRfytzZUCquSY6+PAmW4W0FaSSKNLEIgQbr2\nf2t38m27BLbvmoNjCCSP2FiSgsm6YL6us71T1Oum1QtIuWa8nJJEPbILgQS5ypaONuzdQFVKyatK\nRBpZh0Dyi40lKSp4GZVWlaygOcvllCRYhECCUDV21lldlTQfzedVJYqiHtmIQPKOFSWp7j5v41Up\n8QKSnqP5JOSQwZJEGlmKQILFrK5KOeE/COxFIPlIeElq6DJY41VJCAmVqJKRkkQ9sheBBFmS3ZTB\n26okM4cMIo2sRiB5SnhJSkBnVZJwBZItAcx2O8RHIEGQ9Pesc74qUYlqoB7ZjkDyl7SSlNUdVMOq\npNx61bb330hPSSKNHEAgQYTM7+ed3wye5vm63Q/u3q12KzujCGgIgeQ1aSUpc9Jm8OIfzff04qeD\nD+743h2tC1vzHJQLqEduIJBgXq7HHVm3LzyIorCEHTt4zORoMpLrrB1p5AwCyXfGS5Kew/ekVaVK\nYSUyvn8PMIVAgi8yqUqd7Z0L/8fC4y8fH9s8dvjD4Tf3vzmrddbnV34+zcDKKpGrcipJ1COXEEgw\nWZL0n02eviq99vxrwQdjx49d9PVFidOISpQeaeQYAgnG6E+jQMqqdOe37xw+Nzx+0vg5i+aoKQ3/\ndZ9ziItkURuBBKUErCTpl7gqzVk6J9l39GRqThvqkXsIJJhhqh5FNVqVOts7r/9v1zf6XaKV6KUt\nLzX61x2TVUkijZxEIOHvdJYkCWkUarQqvX/o/eN/OD5x6sR3/++7C/9x4cTPTBztK8sqkeaj+RzW\n/Wg3aeQkAglooCq98/o7g+8M3vzPNyulWv7S8vy/PH/vk/dO/cLU6NeMtkqk52g++dKXpGJPUT2R\n4YggBYGEy7oOd219ZOuqzaty/S6i6lFU1aoU3sw0sHzL8qtuuCr4+Kobrvp026d3/+/d979wv/J7\nt4JOTNY5jEDCxxR7impzjs8vNo0C0RuzVvX8Pz8fzZtJMyYdeeUIuxUalbgkkUZuI5CAcnVjKRBW\noju/fWfirXcAQmNMDwCydB3u6mjryOnJhdejqBpp9PTip4N/1r2y7gvLv6CUGjOW36OGBSWpob9C\nPXIeDQmaWJRGdYWzc0NnhprGNF1767Vmx+MD0sgHvLNDuVxLklMG1bF/P3bb125r+kST6aFYKUFJ\ngttoSNDBpXqklBr448DIxZE/bv/jratuDbaAI1fUI0/QkFBFtiXJxjSqcUHSmu1rzr1/7sK5C3dt\nvuvWb9wa8wnjH83nlTgliTTyBw0JaMyV11155XVXmh4F4CAaEqrLqiTZWI8CG/ZuKOtJ615Zx8VG\nmatdkqhHXiGQkCN70ygUxhJRpB9p5BsCCaNiux30YLsdAgQS8uJAPYJB1CMPEUioJXFJIo3QkLKS\nRBr5iUACAIhAIKGOBCWJelSGo/niCEsS9chbBBIyRhpVGjg4YHoI1iCNfEYgoT622wHQgEBClqhH\nSGPv43vnLZ9nehQwhkBCLHFKEmmENDrbOwt9BdOjgEkEErJBGiETqzav2vrIVtOjgBkEEuJiJQn5\noR5BEUjIBPUIaZSlESXJWwQSGlC1JJFGtXEREhATgQTki4uQaqs6WUdJ8hOBhMaUlSTqEdJg6QhR\nBBKSI42QH0qShwgkNIztdsgE9QhlCCQkRD1CGnHSiJLkGwIJScxbPo+1eiS267FddCNUIpCQULGn\naHoIsFX8Hx5KklcIJDRs6yNbV21eZXoUsBVLRxgNgYTGhGlU6Ct0tneaHo50XBVbJkEaUZL8QSAB\nOWKlDYiPQEIDyibrKEloSOLJOkqSJwgkxMXSEdJg6Qh1EUhIhZKkwdGjR4cvDr/3/nsD7w4MDAz0\n9/ef/fCs6UHpRknywTjTA4AdqEdmnTh+Qil13XXXKaU+PPfhqXdODU8dnnr1VNPjiot6hDhoSKiv\ndhpRkjS7ctKVSqkzZ86YHkhcWaURJcl5BBJgjaNHj4YfN6kmgyMB8kAgoY44k3WUpKpyugjp/N/O\nK6WunHxl5s+ch2wn6yhJbiOQUAtLR2mkvwip/63+ygfPnD7T3Nw8bdq0lE+uAUtHaAiBhGxQkjLU\n/1Z/8I9SasLVE86fOR9+6uRfTzaNafrU7E+NHTPW3ABNoiQ5jF12GBX1SL8ghFo+0xL88djBYxOu\nnhBmUuvCVqXUu++9+9eTf71m1jXjx403ONS6qEdoFA0J1SVII0pSYmElavlMS5hGUUEsvf7a60qp\nGdNnDI8Mv/3Xt0cujWgfaVy5phElyVU0JMCksko0muAKJKVUkEnTZk879+G5s2fPfuKqT+Q9QkAb\nGhKqSDxZR0mKqW4lUkodO3isdWHrddddN755/FtvvTV8cVgp1bqwtXVh6/sn3z9/5vyFCxf0jjou\nJuuQDIGEciwd5acsh+oWI6VUSZWGLwyPjIyMjFyeoJv6yakTrp7w7vF38xxsQnrSiFk7JxFIyBgl\nqapGcyjUpJqumHDF9BnTJ1wxIXhk5NLI0PmhcePHff6Wz7/+2uvBJB7gAAIJH0M9ykpwVWycqbm6\nZl0za2ho6MwHZ4aGhs4NnRs4NdDc3NzyyZYxTWOCGTw5maRzso6S5B4CCZdllUaUJKXUwMGBZJUo\nECwgBR+PHTt21jWzJk2cNDIycmnk0tVTr579qdnjx17e8x1kkvFYYukIKRFI+Du6USai17QmrkRV\nNTc3T548efLkyeHcXZS0qqQHJckxBBJy4VtJSrBbIQ8GqxL1COkRSFCKepSChByKMlKVDKYRJckl\nBBLySiO3S1ImuxVGE11ASkbIqhLQEAIJaIy0SjQabVXJ+GQdJckZBJLvcp2sK/QVdj22K6cn1yzX\nSpSfvKuS8TSCSwgkr2lYOir2FHN9/ryFr+YNVaKcjuZLxocNeJQkN3BzVaC64BU8WMvZrXY39HcH\nDg6oZQm/b/oFpKrCTMrwyalHyBYNyV96dtZZt7UhqESvv/Z6UCxMDydLYVXKpC1JSyNKkgNoSJ5i\nn3elzAuETMG/YJC4iZ9EWhrBDTQk5E54SRJVifrf6tczBif3hVOSbEdD8hH1SH3Uh5SwShQcVa5H\n4qpEPUJOaEjeMZJGokpStA+JSiMjGq1KwtOIkmQ1GhJ8IbMSSZDJqhKQHg3JLwYn60yVpLIlIsmv\nuTlt+I4pTlUSXo8ClCR70ZA84tvSkaldc6Kuim1I7apkRRrBajQk6KOnJBnfNTdwcED/N82QAxvw\nKEmWoiH5wod65MmFRBpUViXqETSgIXlBThrlUZKMV6JMmF1AqiqsSjamESXJRjQk2Ipdcxoku5Uf\nkAwNyX1y6lEgfUmyZdecG4J61Luzt3dnr+mxNIaSZB0akuOkpVEaVCL9wsm6pSuXKqV6d/YGHwB5\nIJBgQKGv0NHWsWHvhphf7/xuBYELSFUtXbk06EnEEvJAILnM9npkYyWy9yKkSlX3MoRVSdkQS8Gs\n3UNPPGR6IIiFQHKW8DSqUZJszKFQmqP5RKm9s44ZPOSBQIIgzk/NOcaKGTxKkkUIJDcJr0eBsCRZ\nXYnSE7iAFP/CI6oSMsS2bwdZkUYhNnBLk+Ay2KAqid0Xzv5vW9CQXGNFGoWvXEFJat1LFFmPqoT0\nCCToE+YQr1khbWeWx5TyLkFiV5VYSbICgeQUsfWoxotUo9ckOUbnmeV1ZXLPOqoSEiOQ3CEwjahE\n3hJYlShJ8hFIyEVDL0aFvkL3o913b7o750HlzuqrYjO/pTdVCY0ikBwhpB4lrkTFnqLalMOA9Gr0\nqlg5G77zO2BCVFWiJAlHILnAeBqln5pzpiShElUJMRFISCXDN79ulCQb6Tl/T0hVoiRJRiBZz0g9\nYreCM3SeBktVQm0Ekt30p1F+b3J92/8tYQFp12O79J9NbrwqUZLEIpAQC5XIScWeotps4PvadYYF\ntCGQLKahHmnOId9Kklk6J+uqMjiDR0mSiUCyVd5pxLvXRtl1EZLxNAoZn8GDHNztGx8T3LM5eNNq\n5AWi0FfobO/U/33TGzg4EP+LJSwgyRH8sGm+WTi3ABeIhmSlPOoR71L9IaceRVGVQCDZJ9s0Erhb\ngZWkXMlMo4DmVSVWkqQhkDwlMIf8Ie3ICWmoSt4ikCyTvh5Z8avudkkyeOSE5HoUpa0qUZJEIZBs\nkiaNqESwJY1CVCXfEEiOszeH3C5JiElDVaIkyUEgWaPResRbS5lMLSBZV4+iqEqeIJDsED+N7K1E\nlSwqSfGvijWygGR1GgVyrUqUJCEIJHfwFtKgRo/mQzJUJbcRSBaoXY9cqkSVLCpJYjlQj6I4w8Jh\nBJI4TU1NwQc7duxYu3ZtmEZnPzhb/M/iLV+6Jfis2znkKv13DHIsjUKZVyVm7SQgkAQJoygUvdfW\nf+z/j+YJzcq/qTlKEqriDAv3cHNV6WbeNFMpdeTPR87+v7Pvn3rf4G1PYR1X61FUhjdm5XarxhFI\nUlTWo4GDA0Ea9e7sPfKfR2Z8csa0a6b5GUX23gLcIB/SKBRkkub7hSNzBJIdlq5cOnH6xIkzJ5oe\nCJLjyIlcZVKVKElmEUhyFXuKpocgiOSSJPBoPq/qURRVyWoEkmh7H9+rlJpy7RTTA0EtDR3Np4G3\naRRIWZUoSQYRSFKUSqWqj+99fO/uB3drHoxMkksSpKEq2YhAskCpVOpo67h0+pLpgSA5bQtIntej\nqMRViZJkCoEkSKlUKutJe/bs2bNnT29v7549e368/sc/XvPjDz74YHBw0NQIjSv0FXY9tsv0KOQi\njSpRlSzChbHiBJm0b9++7du3/+xnPwsfv+KKK2bNmvXUkqe+/oOvT5nv76pSsaeoNpkeBKyS4G5D\n3LjBCAJJqEWLFi1atCj6yAsvvHDDDTe8+eabwRVLvBG2iJ4jJ6hHtXFjVvmYsrPG8PDw6dOn1Ucz\nex1tHR1tHaYHZYCNWxs0HDlBGsXR0KoSK0n6EUgWWLVq1fz58z/44IO+vr4FCxZ85StfUZFYMj06\nwDKsKonVNNpuY9jCzxk8zbdb7WzvXPfKuqqfOvTLQzcuu7H2X897ix31KJk4M3jdj3ZHV5LWL1jP\na2Z+aEjWoyqZVfeq2LwXkEijxDK8MSsyQSA5wrdVJYtWknJdQNr12C7SKKXaM3isJOlEILmDquQh\nbniYCaqSEASSa/ypShaVpJwwWZet0aoSJUkbAslBYVVqmdhieiy+y287A2mUh/hV6S9/+cuaNWvu\nu++++fPnL1u27OWXX9YwPOcRSM4qlUrrF6x3uypRkpCHyqoUlKQdO3YopZqamorF4re+9a1t27b9\n5Cc/+fOf/3zTTTctWbJk27Zt5obsCALJZVQlV1GP8hatSr07e7sf7Y5+9nOf+9yLL744derU4I9b\ntmyZO3fuww8/fPLkSQNjdQiB5D63q5LZkmTkaD7SSJswk2beNFNVbPF/9dVXw48XL148PDz8yiuv\naB6hYwgkL1CVclL7IiTOLHfb4sWLw4/HjRunlLp48aK54biAQPKIq1XJq5Uk6pEpQUkaTbFYVEq1\ntbXpGo6buNu3X4K7nvh5tyEHkEZmjXbV1+HDh3/3u9+tXr16/vz5mofkGALJR2EsOfPqVugraL67\nHTwx5dopux/cXftrvvGNb3zpS1/6/ve/r2dIDiOQ/FUqlahK+cl8AYl6pF/t+e3gjd3GjRubm5t/\n/vOfT548Wde4nEUgec2lquR2SSKNdApzKDrFXdVTTz319ttv//rXvw42NSAl/iOCqgQoVZFDocpY\nCh554YUXDh069IMf/CB4sK+v7/z587fffrum4bqIXXZQypUbs8rZbpftkRPUo7wFt38sfaTq15Qi\nlFIvv/xyX1/fM888E37BgQMHhoaGNI3YUTQkXEZVakiNq2IzPHKCNMrPaJWorgMHDqxevbq9vf2+\n++4LHrl06dIf/vCHo0ePZjxEzxBI+JhwdqLrcFf/UL/p4TRM50rSwMEBtUzD90H2gihKfPbrihUr\n+vv7f/rTn0YfnDJlSnNzcwaD8xiBhCqoSkJQj7KVuBKVOX78eBbDQTkCCdXZW5WMb7fLagGJNMpK\nVjmEvBFIqIWqlECuZ5ajISmn5qAZgYQ6bKxKhb5C96Pdd2+62/RAkqMepUElshSBhFisq0rFnqLa\nZHoQSZFGiVGJrEYgIS67qpKpksSRE0ZQidxAIKExFlWlXEtSfkfzUY/iI4ccQyChYZxhoXK7CIk0\niompOScRSEhI/o1Zje//RuaoRG4jkJCKRTN4GqRcQKIejaZlYsv6BesVOeQ6AglpSa5KFpUk0qgq\nKpFXCCRkg6qUxq7HdvHfLYoc8hPHTyAzMs+wkHMmRQ3FnqLpIUgR5yQIuIqGV6b4sAAABG1JREFU\nhIx5W5USLyAxWaeoRFBKEUjIg7RLaCWvJHmeRuQQoggk5MXhqpTfVbH+4EIiVCKQkCM5VSnbklR5\nVWyyIyc8rEdUItRAICF3DlelUIIjJ3xLIyoR6iKQoIOEqiR5JclhVCLERyBBHx+qUkzO1yNyCAkQ\nSNDKbFXKqSQ1uoDkdhoxNYfECCQY4FhV4sxyRSVCFggkmGHqDAvjK0nu1SMqEbJCIMEkyTdmHU2a\ni5BcSiMqETJHIME8zTN4KUtS2UVIvp1ZTg4hPwQSRLCxKjXK9nrE1BzyRiBBEG1VSf9Kkr1pRCWC\nNgQSZPGhKtmCSgTNCCRIpKEqZVKSYi4g2VWPqEQwhUCCUBLuNpQJW9KIHIJxBBJEy7UqGb8mSQim\n5iAEgQTprK5KkusRlQjSEEiwQ05VqdBX6H60++5Nd8f8+uhVsXUXkMSmEZUIMhFIsEZOVanYU1Sb\n4n5x5dF8FqESQTgCCZax5cascuoROQRbEEiwT7Y3Zk22taH2kRNC0oipOdiFQIKtzF5CK/nICSoR\nLEUgwW6ZzOBlu//bVD1qmdiyfsF6RQ7BWgQSrCfqbkNG0ohKBDcQSHBEyqrUUEkScuQEOQTHEEhw\nR95Vqe7RfNrqEbsV4CQCCa5JXJXqlqTaFyFpSCMqEdxGIMFB+leVdj22K7/vxW4FeIJAgrMSVKU4\nK0lVF5CKPUW1Odkwa6ESwSsEElym58asmU/WkUPwE4EE9zVUlRq9JinbNGK3AnxGIMELws+woBIB\nikCCV2JWpRolqWwBKWU9IoeAKAIJfsmwKqVJI6bmgEoEEnxUtypVlqTgqthjB4+l+b5UIqAGAgme\navQMi7KrYhutR1QioC4CCV6rcQltZUkKF5DipxGVCIiPQAKyP4WWHAISGGN6AIAIpVKpVCqFQRIo\n9BU62zvLvrKzvfPbv/r2c+uee7Pvzcrn6Wjr6GjrKH0kxxEDzqEhAZdVrUphJi2Yt+DJ9icXb1j8\nzP985uShk5cuXgq/hkoEpEcgAR8T3exQpuepHqXUP/zTP0yfPv3koZPBg+xWALJCIAFVhFUpqthT\njP5xx4M7FFEEZKeJXyegUtWGVGbPnj1Lly7VMBjAE2xqABr23HPPmR4C4CACCQAgAoEEABCBQAKq\nqLG2yrIrkBMCCQAgAtu+geoqL0iiGwG5IpCAWgghQBum7ICGDQ4OKqUuXrxoeiCAU7gwFmjAvffe\ne+LEiQMHDoyMjIwfP/7222+fO3fus88+a3pcgAsIJACACEzZAQBEIJAAACIQSAAAEQgkAIAIBBIA\nQAQCCQAgAoEEABCBQAIAiEAgAQBEIJAAACIQSAAAEQgkAIAIBBIAQAQCCQAgAoEEABCBQAIAiEAg\nAQBEIJAAACIQSAAAEQgkAIAIBBIAQAQCCQAgAoEEABCBQAIAiEAgAQBEIJAAACIQSAAAEQgkAIAI\nBBIAQAQCCQAgAoEEABCBQAIAiEAgAQBEIJAAACIQSAAAEQgkAIAIBBIAQIT/DwHyO/eww8neAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[node,elem] = cubemesh([-1,1,-1,1,-1,1],2);\n",
    "bdFlag = setboundary3(node,elem,'Dirichlet','x==1','Neumann','x~=1');\n",
    "figure(2); clf;\n",
    "showmesh3(node,elem);\n",
    "display(elem); display(bdFlag);\n",
    "findnode3(node,[1,2,3,4,5,7,8]);\n",
    "display('change to ascend ordering');\n",
    "[elem,bdFlag] = sortelem3(elem,bdFlag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `bdFlag` to find the boundary nodes, edges and faces. To find the outwords normal direction of the boundary face, we use `gradbasis3` to get `Dlambda(t,:,k)` which is the gradient of $\\lambda_k$. The outward normal direction of the kth face can be obtained by `-Dlambda(t,:,k)` which is independent of the ordering and orientation of `elem`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "DirichletFace =\n",
      "\n",
      "  2�3 uint32 matrix\n",
      "\n",
      "   2   3   7\n",
      "   2   6   7\n",
      "\n",
      "\n",
      "NeumannFace =\n",
      "\n",
      "  10�3 uint32 matrix\n",
      "\n",
      "   1   2   3\n",
      "   1   2   6\n",
      "   1   3   4\n",
      "   1   4   8\n",
      "   1   5   6\n",
      "   1   5   8\n",
      "   3   4   7\n",
      "   4   7   8\n",
      "   5   6   7\n",
      "   5   7   8\n",
      "\n",
      "\n",
      "bdFaceOutDirec =\n",
      "\n",
      "     0     0    -1\n",
      "     0    -1     0\n",
      "     0     0    -1\n",
      "    -1     0     0\n",
      "     0    -1     0\n",
      "    -1     0     0\n",
      "     0     1     0\n",
      "     0     1     0\n",
      "     0     0     1\n",
      "     0     0     1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Dlambda = gradbasis3(node,elem);\n",
    "T = auxstructure3(elem);\n",
    "elem2face = T.elem2face; \n",
    "face = T.face;\n",
    "NF = size(face,1);\n",
    "if ~isempty(bdFlag)\n",
    "    % Find Dirchelt boundary faces and nodes\n",
    "    isBdFace = false(NF,1);\n",
    "    isBdFace(elem2face(bdFlag(:,1) == 1,1)) = true;\n",
    "    isBdFace(elem2face(bdFlag(:,2) == 1,2)) = true;\n",
    "    isBdFace(elem2face(bdFlag(:,3) == 1,3)) = true; \n",
    "    isBdFace(elem2face(bdFlag(:,4) == 1,4)) = true;\n",
    "    DirichletFace = face(isBdFace,:);\n",
    "    % Find outwards normal direction of Neumann boundary faces\n",
    "    bdFaceOutDirec = zeros(NF,3);\n",
    "    bdFaceOutDirec(elem2face(bdFlag(:,1) == 2,1),:) = -Dlambda(bdFlag(:,1) == 2,:,1);\n",
    "    bdFaceOutDirec(elem2face(bdFlag(:,2) == 2,2),:) = -Dlambda(bdFlag(:,2) == 2,:,2);\n",
    "    bdFaceOutDirec(elem2face(bdFlag(:,3) == 2,3),:) = -Dlambda(bdFlag(:,3) == 2,:,3);\n",
    "    bdFaceOutDirec(elem2face(bdFlag(:,4) == 2,4),:) = -Dlambda(bdFlag(:,4) == 2,:,4);\n",
    "end\n",
    "% normalize the boundary face outwards direction\n",
    "vl = sqrt(dot(bdFaceOutDirec,bdFaceOutDirec,2));\n",
    "idx = (vl==0);\n",
    "NeumannFace = face(~idx,:);\n",
    "bdFaceOutDirec(idx,:) = [];\n",
    "vl(idx) = [];\n",
    "bdFaceOutDirec = bdFaceOutDirec./[vl vl vl];\n",
    "display(DirichletFace);\n",
    "display(NeumannFace);\n",
    "display(bdFaceOutDirec);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "308px",
    "width": "416px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_position": {
    "height": "359px",
    "left": "0px",
    "right": "1040px",
    "top": "106px",
    "width": "240px"
   },
   "toc_section_display": "none",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
